X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FEEOver.fs;h=489105940e91a575522666029c2c4dc7abf57d36;hp=e366b72889b6d8b1a95d056b0fe425cb9c2c5d53;hb=dec96d50e56e1932bbfa09e6bedf90d6b707ccbd;hpb=ba64921fb9a0c36cd8cf802cbf1b2c0f79bc25f6 diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index e366b72..4891059 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -55,6 +55,16 @@ let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (b let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff +#if DEBUG_LOG + printf "\t\t--- debug istanpt with (x,y)=(%f, %f), A1=%f, B1=%f\n" x y a1 b1 + printf "theta=%f\n" theta + printf "eps_Radian=%f\n" eps_radian + printf "(x1, y1)=(%f, %f)\n" x1 y1 + printf "(x2, y2)=(%f, %f)\n" x2 y2 + printf "test1=%f\n" test1 + printf "test2=%f\n" test2 +#endif + if test1 * test2 > 0.0 then TANGENT_POINT else INTERSECTION_POINT @@ -102,6 +112,9 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: if area1 < 0.0 then +#if DEBUG_LOG + printf "TWO area1=%f\n" area1 +#endif area1 <- area1 + a1 * b1 let cosphi = cos (phi_1 - phi_2) @@ -161,6 +174,9 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: let mutable area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr)) if area2 < 0.0 then +#if DEBUG_LOG + printf "TWO area2=%f\n" area2 +#endif area2 <- area2 + a2 * b2 area1 + area2 @@ -172,9 +188,12 @@ let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) for i in 0..2 do if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT then - tanpts <- tanpts + 1; - tanindex <- i; - + tanpts <- tanpts + 1 + tanindex <- i +#if DEBUG_LOG + printf "tanindex=%d\n" tanindex +#endif + if tanpts <> 1 then -1.0 @@ -205,25 +224,40 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( xint.[i] <- if xint.[i] < 0.0 then -a1 else a1 theta.[i] <- if yint.[i] < 0.0 then 2.0 * Math.PI - acos (xint.[i] / a1) else acos (xint.[i] / a1) +#if DEBUG_LOG + for k in 0..3 do + printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k] +#endif + for j in 1..3 do let tmp0 = theta.[j] let tmp1 = xint.[j] let tmp2 = yint.[j] let mutable k = j - 1 + let mutable k2 = 0 while k >= 0 do if theta.[k] <= tmp0 then + k2 <- k + 1 k <- -1 else theta.[k+1] <- theta.[k] xint.[k+1] <- xint.[k] yint.[k+1] <- yint.[k] - k <- k - 1 + k <- k - 1 + k2 <- k + 1 + + theta.[k2] <- tmp0 + xint.[k2] <- tmp1 + yint.[k2] <- tmp2 - theta.[k+1] <- tmp0 - xint.[k+1] <- tmp1 - yint.[k+1] <- tmp2 + +#if DEBUG_LOG + printf "AFTER sorting\n" + for k in 0..3 do + printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k] +#endif let area1 = 0.5 * abs ((xint.[2] - xint.[0]) * (yint.[3] - yint.[1]) - (xint.[3] - xint.[1]) * (yint.[2] - yint.[0])) @@ -267,20 +301,36 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( if area5 < 0.0 then +#if DEBUG_LOG + printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2 +#endif area5 <- area5 + area_2 if area4 < 0.0 then +#if DEBUG_LOG + printf "\n\t\t-------------> area4 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area4 area_2 +#endif area4 <- area4 + area_2 if area3 < 0.0 then +#if DEBUG_LOG + printf "\n\t\t-------------> area3 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1 +#endif area3 <- area3 + area_1 if area2 < 0.0 - then + then +#if DEBUG_LOG + printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1 +#endif area2 <- area2 + area_1 +#if DEBUG_LOG + printf "\narea1=%f, area2=%f area3=%f, area4=%f, area5=%f\n\n" area1 area2 area3 area4 area5 +#endif + area1 + area2 + area3 + area4 + area5 @@ -463,11 +513,9 @@ let private biquadroots (p: float[]) (r: float[,]) = quad () -open Types - -let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = - let { cx = h1; cy = k1; a = a1; b = b1; alpha = phi_1 } = e1 - let { cx = h2; cy = k2; a = a2; b = b2; alpha = phi_2 } = e2 +let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = + let h1, k1, a1, b1, phi_1 = e1.Cx, e1.Cy, e1.A, e1.B, e1.Alpha + let h2, k2, a2, b2, phi_2 = e2.Cx, e2.Cy, e2.A, e2.B, e2.Alpha if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS then @@ -480,6 +528,10 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = let sinphi = sin phi_1 (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI) +#if DEBUG_LOG + printf "H2_TR=%f, K2_TR=%f, PHI_2R=%f\n" h2_tr k2_tr phi_2r +#endif + let cosphi = cos phi_2r let cosphi2 = cosphi ** 2.0 let sinphi = sin phi_2r @@ -507,6 +559,11 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc) |] +#if DEBUG_LOG + for i in 0..4 do + printf "cy[%d]=%f\n" i cy.[i] +#endif + let py = Array.zeroCreate 5 let r = Array2D.zeroCreate 3 5 @@ -516,6 +573,10 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = for i in 0 .. 3 do py.[4-i] <- cy.[i] / cy.[4] py.[0] <- 1.0 +#if DEBUG_LOG + for i in 0..4 do + printf "py[%d]=%f\n" i py.[i] +#endif biquadroots py r 4 @@ -544,14 +605,27 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = else 0 +#if DEBUG_LOG + printf "nroots = %d\n" nroots +#endif + let ychk = [| for i in 1 .. nroots do if abs r.[2, i] < EPS then yield r.[1, i] * b1 +#if DEBUG_LOG + printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1 +#endif |] Array.sortInPlace ychk +#if DEBUG_LOG + printf "nychk=%d\n" ychk.Length + for j in 0 .. ychk.Length - 1 do + printf "\t j=%d, ychk=%f\n" j ychk.[j] +#endif + let nychk = Array.length ychk let mutable nintpts = 0 @@ -562,30 +636,61 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = let mutable i = 0 while returnValue = 0.0 && i < nychk do +#if DEBUG_LOG + printf "------------->i=%d (nychk=%d)\n" i nychk +#endif + if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0) then +#if DEBUG_LOG + printf "check intersecting points. nintps is %d" nintpts +#endif + let x1 = if abs ychk.[i] > b1 then 0.0 else a1 * sqrt (1.0 - (ychk.[i] * ychk.[i]) / (b1 * b1)) let x2 = -x1 +#if DEBUG_LOG + printf "\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) + 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) +#endif + if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS then nintpts <- nintpts + 1 +#if DEBUG_LOG + printf "first if x1. acc nintps=%d\n" nintpts +#endif if nintpts > 4 then returnValue <- -1.0 else xint.[nintpts-1] <- x1 yint.[nintpts-1] <- ychk.[i] +#if DEBUG_LOG + printf "nintpts=%d, xint=%f, x2=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i] +#endif if returnValue <> -1.0 && abs (ellipse2tr x2 ychk.[i] aa bb cc dd ee ff) < EPS && abs (x2 - x1) > EPS then nintpts <- nintpts + 1 +#if DEBUG_LOG + printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS +#endif if nintpts > 4 then returnValue <- -1.0 else xint.[nintpts-1] <- x2 yint.[nintpts-1] <- ychk.[i] + +#if DEBUG_LOG + printf "nintpts=%d, x1=%f, xint=%f, i=%d, yint=%f\n" nintpts x1 x2 i ychk.[i] +#endif + +#if DEBUG_LOG + else + printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1] +#endif i <- i + 1 @@ -596,9 +701,17 @@ let EEOverlapArea (e1: Ellipse) (e2: Ellipse) : float = match nintpts with | 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 | 2 -> match istanpt xint.[0] yint.[0] a1 b1 aa bb cc dd ee ff with - | 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 - | INTERSECTION_POINT -> twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff + | TANGENT_POINT -> +#if DEBUG_LOG + printf "one point is tangent\n" +#endif + nointpts a1 b1 a2 b2 h1 k1 h2 k2 phi_1 phi_2 h2_tr k2_tr aa bb cc dd ee ff + + | INTERSECTION_POINT -> +#if DEBUG_LOG + printf "check twointpts\n" +#endif + twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff | 3 -> threeintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff | 4 -> fourintpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff | _ -> -1.0 -