Add report.
[master-thesis.git] / Parasitemia / ParasitemiaCore / EEOver.fs
1 module ParasitemiaCore.EEOver
2
3 open System
4
5 let private EPS = 1.0e-7
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 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) =
73 if abs x.[0] > a1
74 then x.[0] <- if x.[0] < 0.0 then -a1 else a1
75
76 let mutable theta1 =
77 if y.[0] < 0.0
78 then 2.0 * Math.PI - acos (x.[0] / a1)
79 else acos (x.[0] / a1)
80
81 if abs x.[1] > a1
82 then x.[1] <- if x.[1] < 0.0 then -a1 else a1
83
84 let mutable theta2 =
85 if y.[1] < 0.0
86 then 2.0 * Math.PI - acos (x.[1] / a1)
87 else acos (x.[1] / a1)
88
89 if theta1 > theta2
90 then
91 let tmp = theta1
92 theta1 <- theta2
93 theta2 <- tmp
94
95 let xmid = a1 * cos ((theta1 + theta2) / 2.0)
96 let ymid = b1 * sin ((theta1 + theta2) / 2.0)
97
98 if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0
99 then
100 let tmp = theta1
101 theta1 <- theta2
102 theta2 <- tmp
103
104 if theta1 > theta2
105 then
106 theta1 <- theta1 - 2.0 * Math.PI
107
108 let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0
109
110 let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0]))
111
112 if area1 < 0.0
113 then
114 #if DEBUG_LOG
115 printf "TWO area1=%f\n" area1
116 #endif
117 area1 <- area1 + a1 * b1
118
119 let cosphi = cos (phi_1 - phi_2)
120 let sinphi = sin (phi_1 - phi_2)
121
122 let mutable x1_tr = (x.[0] - h2_tr) * cosphi + (y.[0] - k2_tr) * -sinphi
123 let mutable y1_tr = (x.[0] - h2_tr) * sinphi + (y.[0] - k2_tr) * cosphi
124 let mutable x2_tr = (x.[1] - h2_tr) * cosphi + (y.[1] - k2_tr) * -sinphi
125 let mutable y2_tr = (x.[1] - h2_tr) * sinphi + (y.[1] - k2_tr) * cosphi
126
127 if abs x1_tr > a2
128 then
129 x1_tr <- if x1_tr < 0.0 then -a2 else a2
130
131 if y1_tr < 0.0
132 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
138 then
139 x2_tr <- if x2_tr < 0.0 then -a2 else a2
140
141 if y2_tr < 0.0
142 then
143 theta2 <- 2.0 * Math.PI - acos (x2_tr / a2)
144 else
145 theta2 <- acos (x2_tr / a2)
146
147 if theta1 > theta2
148 then
149 let tmp = theta1
150 theta1 <- theta2
151 theta2 <- tmp
152
153 let xmid = a2 * cos ((theta1 + theta2) / 2.0)
154 let ymid = b2 * sin ((theta1 + theta2) / 2.0)
155
156 let cosphi = cos (phi_2 - phi_1)
157 let sinphi = sin (phi_2 - phi_1)
158 let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr
159 let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr
160
161 if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0
162 then
163 let tmp = theta1
164 theta1 <- theta2
165 theta2 <- tmp
166
167 if theta1 > theta2
168 then
169 theta1 <- theta1 - 2.0 * Math.PI
170
171 let trsign = if theta2 - theta1 > Math.PI then 1.0 else -1.0
172
173 let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr))
174 if area2 < 0.0
175 then
176 #if DEBUG_LOG
177 printf "TWO area2=%f\n" area2
178 #endif
179 area2 <- area2 + a2 * b2
180
181 area1 + area2
182
183 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 =
184 let mutable tanpts = 0
185 let mutable tanindex = 0
186 for i in 0..2 do
187 if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT
188 then
189 tanpts <- tanpts + 1
190 tanindex <- i
191 #if DEBUG_LOG
192 printf "tanindex=%d\n" tanindex
193 #endif
194
195 if tanpts <> 1
196 then
197 -1.0
198 else
199 match tanindex with
200 | 0 ->
201 xint.[0] <- xint.[2]
202 yint.[0] <- yint.[2]
203 | 1 ->
204 xint.[1] <- xint.[2]
205 yint.[1] <- yint.[2]
206 | _ ->
207 ()
208 twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
209
210 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 =
211 let a1b1 = a1 * b1
212 let a2b2 = a2 * b2
213 let area_1 = Math.PI * a1b1
214 let area_2 = Math.PI * a2b2
215
216 let theta = Array.zeroCreate 4
217
218 for i in 0 .. 3 do
219 if abs xint.[i] > a1
220 then
221 xint.[i] <- if xint.[i] < 0.0 then -a1 else a1
222 theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1)
223
224 #if DEBUG_LOG
225 for k in 0..3 do
226 printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
227 #endif
228
229 for j in 1 .. 3 do
230 let tmp0 = theta.[j]
231 let tmp1 = xint.[j]
232 let tmp2 = yint.[j]
233
234 let mutable k = j - 1
235 let mutable k2 = 0
236 while k >= 0 do
237 if theta.[k] <= tmp0
238 then
239 k2 <- k + 1
240 k <- -1
241 else
242 theta.[k+1] <- theta.[k]
243 xint.[k+1] <- xint.[k]
244 yint.[k+1] <- yint.[k]
245 k <- k - 1
246 k2 <- k + 1
247
248 theta.[k2] <- tmp0
249 xint.[k2] <- tmp1
250 yint.[k2] <- tmp2
251
252
253 #if DEBUG_LOG
254 printf "AFTER sorting\n"
255 for k in 0..3 do
256 printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k]
257 #endif
258
259 let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0]))
260
261 let cosphi = cos (phi_1 - phi_2)
262 let sinphi = sin (phi_1 - phi_2)
263
264 let theta_tr = Array.zeroCreate 4
265 let xint_tr = Array.zeroCreate 4
266 let yint_tr = Array.zeroCreate 4
267
268 for i in 0..3 do
269 xint_tr.[i] <- (xint.[i] - h2_tr) * cosphi + (yint.[i] - k2_tr) * -sinphi
270 yint_tr.[i] <- (xint.[i] - h2_tr) * sinphi + (yint.[i] - k2_tr) * cosphi
271
272 if abs xint_tr.[i] > a2
273 then
274 xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2
275
276 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)
277
278 let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0)
279 let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0)
280
281 let mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0
282
283 if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0
284 then
285 area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs (xint.[0] * yint.[1] - xint.[1] * yint.[0]))
286 area3 <- 0.5 * (a1b1 * (theta.[3] - theta.[2]) - abs (xint.[2] * yint.[3] - xint.[3] * yint.[2]))
287 area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1]))
288
289 if theta_tr.[3] > theta_tr.[0]
290 then
291 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]))
292 else
293 area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3]))
294 else
295 area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1]))
296 area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3]))
297 area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0]))
298 area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2]))
299
300 if area5 < 0.0
301 then
302 #if DEBUG_LOG
303 printf "\n\t\t-------------> area5 is negative (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2
304 #endif
305 area5 <- area5 + area_2
306
307 if area4 < 0.0
308 then
309 #if DEBUG_LOG
310 printf "\n\t\t-------------> area4 is negative (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2
311 #endif
312 area4 <- area4 + area_2
313
314 if area3 < 0.0
315 then
316 #if DEBUG_LOG
317 printf "\n\t\t-------------> area3 is negative (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1
318 #endif
319 area3 <- area3 + area_1
320
321 if area2 < 0.0
322 then
323 #if DEBUG_LOG
324 printf "\n\t\t-------------> area2 is negative (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1
325 #endif
326 area2 <- area2 + area_1
327
328 #if DEBUG_LOG
329 printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5
330 #endif
331
332 area1 + area2 + area3 + area4 + area5
333
334 let private quadroots (p: float[]) (r: float[,]) =
335 let mutable b = -p.[1] / (2.0 * p.[0])
336 let c = p.[2] / p.[0]
337 let mutable d = b * b - c
338
339 if d >= 0.0
340 then
341 if b > 0.0
342 then
343 b <- sqrt d + b
344 r.[1, 2] <- b
345 else
346 b <- -sqrt d + b
347 r.[1, 2] <- b
348 r.[1, 1] <- c / b
349 r.[2, 1] <- 0.0
350 r.[2, 2] <- 0.0
351 else
352 d <- sqrt -d
353 r.[2, 1] <- d
354 r.[2, 2] <- -d
355 r.[1, 1] <- b
356 r.[1, 2] <- b
357
358 let private cubicroots (p: float[]) (r: float[,]) =
359 if p.[0] <> 1.0 then
360 for k in 1..3 do
361 p.[k] <- p.[k] / p.[0]
362 p.[0] <- 1.0
363 let s = p.[1] / 3.0
364 let mutable t = s * p.[1]
365 let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3])
366 t <- (t - p.[2]) / 3.0
367 let mutable c = t * t * t
368 let mutable d = b * b - c
369
370 if d >= 0.0
371 then
372 d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0)
373 if d <> 0.0
374 then
375 if b > 0.0
376 then b <- -d
377 else b <- d
378 c <- t / b
379 d <- sqrt(0.75) * (b - c)
380 r.[2, 2] <- d
381 b <- b + c
382 c <- -0.5 * b - s
383 r.[1, 2] <- c
384 if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0
385 then
386 r.[1, 1] <- c
387 r.[2, 1] <- -d
388 r.[1, 3] <- b - s
389 r.[2, 3] <- 0.0
390 else
391 r.[1, 1] <- b - s
392 r.[2, 1] <- 0.0
393 r.[1, 3] <- c
394 r.[2, 3] <- -d
395 else
396 if b = 0.0
397 then d <- (atan 1.0) / 1.5
398 else d <- atan ((sqrt -d) / (abs b)) / 3.0
399
400 if b < 0.0
401 then b <- 2.0 * (sqrt t)
402 else b <- -2.0 * (sqrt t)
403
404 c <- (cos d) * b
405 t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c
406 d <- -t - c - s
407 c <- c - s
408 t <- t - s
409
410 if abs c > abs t
411 then
412 r.[1, 3] <- c
413 else
414 r.[1, 3] <- t
415 t <- c
416
417 if abs d > abs t
418 then
419 r.[1, 2] <- d
420 else
421 r.[1, 2] <- t
422 t <- d
423
424 r.[1, 1] <- t
425 for k in 1..3 do
426 r.[2, k] <- 0.0
427
428 let private biquadroots (p: float[]) (r: float[,]) =
429 if p.[0] <> 1.0
430 then
431 for k in 1..4 do
432 p.[k] <- p.[k] / p.[0]
433 p.[0] <- 1.0
434 let e = 0.25 * p.[1]
435 let mutable b = 2.0 * e
436 let mutable c = b ** 2.0
437 let mutable d = 0.75 * c
438 b <- p.[3] + b *(c - p.[2])
439 let mutable a = p.[2] - d
440 c <- p.[4] + e * (e * a - p.[3])
441 a <- a - d
442
443 let mutable quadExecuted = false
444 let inline quad () =
445 if not quadExecuted
446 then
447 p.[2] <- c / b
448 quadroots p r
449 for k in 1..2 do
450 for j in 1..2 do
451 r.[j, k+2] <- r.[j, k]
452 p.[1] <- -p.[1]
453 p.[2] <- b
454 quadroots p r
455 for k in 1..4 do
456 r.[1,k] <- r.[1,k] - e
457 quadExecuted <- true
458
459 p.[1] <- 0.5 * a
460 p.[2] <- (p.[1] * p.[1] - c) * 0.25
461 p.[3] <- b * b / -64.0
462 if p.[3] < 0.0
463 then
464 cubicroots p r
465 let mutable k = 1
466 while k < 4 do
467 if r.[2, k] = 0.0 && r.[1, k] > 0.0
468 then
469 d <- r.[1, k] * 4.0
470 a <- a + d
471 if a >= 0.0 && b >= 0.0
472 then
473 p.[1] <- sqrt d
474 elif a <= 0.0 && b <= 0.0
475 then
476 p.[1] <- sqrt d
477 else
478 p.[1] <- -(sqrt d)
479 b <- 0.5 * (a + b / p.[1])
480 quad ()
481 k <- 4
482 k <- k + 1
483
484 if not quadExecuted && p.[2] < 0.0
485 then
486 b <- sqrt c
487 d <- b + b - a
488 p.[1] <- 0.0
489 if d > 0.0
490 then
491 p.[1] <- sqrt d
492 elif not quadExecuted
493 then
494 if p.[1] > 0.0
495 then
496 b <- (sqrt p.[2]) * 2.0 + p.[1]
497 else
498 b <- -(sqrt p.[2]) * 2.0 + p.[1]
499
500 if b <> 0.0
501 then
502 p.[1] <- 0.0
503 else
504 for k in 1..4 do
505 r.[1, k] <- -e
506 r.[2, k] <- 0.0
507 quadExecuted <- true
508
509 quad ()
510
511 // Return a tuple (area, x intersections, y intersections)
512 let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
513 let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha
514 let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha
515
516 if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS
517 then
518 None
519 else
520 let phi_1 = phi_1 % Math.PI
521 let phi_2 = phi_2 % Math.PI
522 let h2_tr, k2_tr, phi_2r =
523 let cosphi = cos phi_1
524 let sinphi = sin phi_1
525 (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI)
526
527 #if DEBUG_LOG
528 printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r
529 #endif
530
531 let cosphi = cos phi_2r
532 let cosphi2 = cosphi ** 2.0
533 let sinphi = sin phi_2r
534 let sinphi2 = sinphi ** 2.0
535 let cosphisinphi = 2.0 * cosphi * sinphi
536 let a22 = a2 ** 2.0
537 let b22 = b2 ** 2.0
538 let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22
539 let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22
540 let tmp2 = cosphi * h2_tr + sinphi * k2_tr
541 let tmp3 = sinphi * h2_tr - cosphi * k2_tr
542
543 let aa = cosphi2 / a22 + sinphi2 / b22
544 let bb = cosphisinphi / a22 - cosphisinphi / b22
545 let cc = sinphi2 / a22 + cosphi2 / b22
546 let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1
547 let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1
548 let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0
549
550 let cy = [|
551 (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff)
552 2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff)
553 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)
554 2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee))
555 a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc)
556 |]
557
558 #if DEBUG_LOG
559 for i in 0..4 do
560 printf "cy[%d]=%f\n" i cy.[i]
561 #endif
562
563 let py = Array.zeroCreate<float> 5
564 let r = Array2D.zeroCreate<float> 3 5
565
566 let nroots =
567 if abs cy.[4] > EPS
568 then
569 for i in 0 .. 3 do
570 py.[4-i] <- cy.[i] / cy.[4]
571 py.[0] <- 1.0
572 #if DEBUG_LOG
573 for i in 0..4 do
574 printf "py[%d]=%f\n" i py.[i]
575 #endif
576 biquadroots py r
577 4
578
579 elif abs cy.[3] > EPS
580 then
581 for i in 0..2 do
582 py.[3 - i] <- cy.[i] / cy.[3]
583 py.[0] <- 1.0
584 cubicroots py r
585 3
586
587 elif abs cy.[2] > EPS
588 then
589 for i in 0..1 do
590 py.[2-i] <- cy.[i] / cy.[2]
591 py.[0] <- 1.0
592 quadroots py r
593 2
594
595 elif abs cy.[1] > EPS
596 then
597 r.[1, 1] <- -cy.[0] / cy.[1]
598 r.[2, 1] <- 0.0
599 1
600
601 else
602 0
603
604 #if DEBUG_LOG
605 printf "nroots = %d\n" nroots
606 #endif
607
608 let ychk = Array.init nroots (fun _ -> Double.MaxValue)
609 let mutable nychk = 0
610 for i in 1 .. nroots do
611 if abs r.[2, i] < EPS
612 then
613 ychk.[nychk] <- r.[1, i] * b1
614 nychk <- nychk + 1
615 #if DEBUG_LOG
616 printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1
617 #endif
618 Array.sortInPlace ychk
619
620 #if DEBUG_LOG
621 printf "nychk=%d\n" ychk.Length
622 for j in 0 .. ychk.Length - 1 do
623 printf "\t j=%d, ychk=%f\n" j ychk.[j]
624 #endif
625
626 let mutable nintpts = 0
627
628 let xint = Array.zeroCreate 4
629 let yint = Array.zeroCreate 4
630
631 let mutable returnValue = 0.0
632
633 let mutable i = 0
634 while returnValue = 0.0 && i < nychk do
635 #if DEBUG_LOG
636 printf "------------->i=%d (nychk=%d)\n" i nychk
637 #endif
638
639 if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0)
640 then
641 #if DEBUG_LOG
642 printf "check intersecting points. nintps is %d" nintpts
643 #endif
644
645 let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1))
646 let x2 = -x1
647
648 #if DEBUG_LOG
649 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)
650 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)
651 #endif
652
653 if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS
654 then
655 nintpts <- nintpts + 1
656 #if DEBUG_LOG
657 printf "first if x1. acc nintps=%d\n" nintpts
658 #endif
659 if nintpts > 4
660 then
661 returnValue <- -1.0
662 else
663 xint.[nintpts-1] <- x1
664 yint.[nintpts-1] <- ychk.[i]
665 #if DEBUG_LOG
666 printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
667 #endif
668
669 if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS
670 then
671 nintpts <- nintpts + 1
672 #if DEBUG_LOG
673 printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS
674 #endif
675 if nintpts > 4
676 then
677 returnValue <- -1.0
678 else
679 xint.[nintpts-1] <- x2
680 yint.[nintpts-1] <- ychk.[i]
681
682 #if DEBUG_LOG
683 printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i]
684 #endif
685
686 #if DEBUG_LOG
687 else
688 printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1]
689 #endif
690 i <- i + 1
691
692 if returnValue = -1.0
693 then
694 None
695 else
696 let area =
697 match nintpts with
698 | 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
699 | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with
700 | TANGENT_POINT ->
701 #if DEBUG_LOG
702 printf "one point is tangent\n"
703 #endif
704 nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff
705
706 | INTERSECTION_POINT ->
707 #if DEBUG_LOG
708 printf "check twointpts\n"
709 #endif
710 twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
711 | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
712 | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff
713 | _ -> -1.0
714
715 if area = -1.0
716 then
717 None
718 elif nintpts = 0
719 then
720 Some (float32 area, [||], [||])
721 else
722 let xTransform : float32[] = Array.zeroCreate nintpts
723 let yTransform : float32[] = Array.zeroCreate nintpts
724 for i in 0 .. (nintpts - 1) do
725 xTransform.[i] <- float32 <| cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1
726 yTransform.[i] <- float32 <| sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1
727 Some (float32 area, xTransform, yTransform)