X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=blobdiff_plain;f=Parasitemia%2FParasitemia%2FEEOver.fs;h=5a3c9e4c580be0c6f6feb74a71c8503a261d1bb1;hp=489105940e91a575522666029c2c4dc7abf57d36;hb=044b0ae69df3ac565432545b2fa934589016f9bd;hpb=dec96d50e56e1932bbfa09e6bedf90d6b707ccbd diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs index 4891059..5a3c9e4 100644 --- a/Parasitemia/Parasitemia/EEOver.fs +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -7,7 +7,7 @@ let private EPS = 1.0e-5 let inline private ellipse2tr (x: float) (y: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : float = aa * x * x + bb * x * y + cc * y * y + dd * x + ee * y + ff -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) = +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) = let a1b1 = a1 * b1 let a2b2 = a2 * b2 let area_1 = Math.PI * a1b1 @@ -34,13 +34,13 @@ let private nointpts (a1: float) (b1: float) (a2: float) (b2: float) (h1: float) type private PointType = TANGENT_POINT | INTERSECTION_POINT let private istanpt (x: float) (y: float) (a1: float) (b1: float) (aa: float) (bb: float) (cc: float) (dd: float) (ee: float) (ff: float) : PointType = - let x = + let x = if abs x > a1 then if x < 0.0 then -a1 else a1 else x - let theta = + let theta = if y < 0.0 then 2.0 * Math.PI - acos (x / a1) else acos (x / a1) @@ -105,7 +105,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: if theta1 > theta2 then theta1 <- theta1 - 2.0 * Math.PI - + let trsign = if (theta2 - theta1) > Math.PI then 1.0 else -1.0 let mutable area1 = 0.5 * (a1 * b1 * (theta2 - theta1) + trsign * abs (x.[0] * y.[1] - x.[1] * y.[0])) @@ -182,7 +182,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: area1 + area2 -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 = +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 = let mutable tanpts = 0 let mutable tanindex = 0 for i in 0..2 do @@ -205,7 +205,7 @@ let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) | 1 -> xint.[1] <- xint.[2] yint.[1] <- yint.[2] - | _ -> + | _ -> () twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff @@ -244,8 +244,8 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( else theta.[k+1] <- theta.[k] xint.[k+1] <- xint.[k] - yint.[k+1] <- yint.[k] - k <- k - 1 + yint.[k+1] <- yint.[k] + k <- k - 1 k2 <- k + 1 theta.[k2] <- tmp0 @@ -258,7 +258,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( 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])) let cosphi = cos (phi_1 - phi_2) @@ -275,7 +275,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( if abs xint_tr.[i] > a2 then xint_tr.[i] <- if xint_tr.[i] < 0.0 then -a2 else a2 - + 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) let xmid = a1 * cos ((theta.[0] + theta.[1]) / 2.0) @@ -286,6 +286,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( if ellipse2tr xmid ymid aa bb cc dd ee ff < 0.0 then area2 <- 0.5 * (a1b1 * (theta.[1] - theta.[0]) - abs (xint.[0] * yint.[1] - xint.[1] * yint.[0])) + area3 <- 0.5 * (a1b1 * (theta.[3] - theta.[2]) - abs (xint.[2] * yint.[3] - xint.[3] * yint.[2])) area4 <- 0.5 * (a2b2 * (theta_tr.[2] - theta_tr.[1]) - abs (xint_tr.[1] * yint_tr.[2] - xint_tr.[2] * yint_tr.[1])) if theta_tr.[3] > theta_tr.[0] @@ -298,9 +299,9 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( area3 <- 0.5 * (a1b1 * (theta.[0] - (theta.[3] - 2.0 * Math.PI)) - abs (xint.[3] * yint.[0] - xint.[0] * yint.[3])) area4 <- 0.5 * (a2b2 * (theta_tr.[1] - theta_tr.[0]) - abs (xint_tr.[0] * yint_tr.[1] - xint_tr.[1] * yint_tr.[0])) area5 <- 0.5 * (a2b2 * (theta_tr.[3] - theta_tr.[2]) - abs (xint_tr.[2] * yint_tr.[3] - xint_tr.[3] * yint_tr.[2])) - + if area5 < 0.0 - then + then #if DEBUG_LOG printf "\n\t\t-------------> area5 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area5 area_2 #endif @@ -321,7 +322,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( 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 @@ -334,7 +335,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( area1 + area2 + area3 + area4 + area5 -let private quadroots (p: float[]) (r: float[,]) = +let private quadroots (p: float[]) (r: float[,]) = let mutable b = -p.[1] / (2.0 * p.[0]) let c = p.[2] / p.[0] let mutable d = b * b - c @@ -342,7 +343,7 @@ let private quadroots (p: float[]) (r: float[,]) = if d >= 0.0 then if b > 0.0 - then + then b <- sqrt d + b r.[1, 2] <- b else @@ -370,7 +371,7 @@ let private cubicroots (p: float[]) (r: float[,]) = let mutable c = t * t * t let mutable d = b * b - c - if d >= 0.0 + if d >= 0.0 then d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0) if d <> 0.0 @@ -427,7 +428,7 @@ let private cubicroots (p: float[]) (r: float[,]) = r.[1, 1] <- t for k in 1..3 do r.[2, k] <- 0.0 - + let private biquadroots (p: float[]) (r: float[,]) = if p.[0] <> 1.0 @@ -436,33 +437,33 @@ let private biquadroots (p: float[]) (r: float[,]) = p.[k] <- p.[k] / p.[0] p.[0] <- 1.0 let e = 0.25 * p.[1] - let b = ref (2.0 * e) - let c = ref (!b ** 2.0) - let mutable d = 0.75 * !c - b := p.[3] + !b *(!c - p.[2]) + let mutable b = 2.0 * e + let mutable c = b ** 2.0 + let mutable d = 0.75 * c + b <- p.[3] + b *(c - p.[2]) let mutable a = p.[2] - d - c := p.[4] + e * (e * a - p.[3]) + c <- p.[4] + e * (e * a - p.[3]) a <- a - d - let quadExecuted = ref false - let quad () = - if not !quadExecuted + let mutable quadExecuted = false + let quad () = + if not quadExecuted then - p.[2] <- !c / !b + p.[2] <- c / b quadroots p r for k in 1..2 do for j in 1..2 do r.[j, k+2] <- r.[j, k] p.[1] <- -p.[1] - p.[2] <- !b + p.[2] <- b quadroots p r for k in 1..4 do r.[1,k] <- r.[1,k] - e - quadExecuted := true + quadExecuted <- true p.[1] <- 0.5 * a - p.[2] <- (p.[1] * p.[1] - !c) * 0.25 - p.[3] <- !b * !b / -64.0 + p.[2] <- (p.[1] * p.[1] - c) * 0.25 + p.[3] <- b * b / -64.0 if p.[3] < 0.0 then cubicroots p r @@ -472,62 +473,62 @@ let private biquadroots (p: float[]) (r: float[,]) = then d <- r.[1, k] * 4.0 a <- a + d - if a >= 0.0 && !b >= 0.0 + if a >= 0.0 && b >= 0.0 then p.[1] <- sqrt d - elif a <= 0.0 && !b <= 0.0 + elif a <= 0.0 && b <= 0.0 then p.[1] <- sqrt d else p.[1] <- -(sqrt d) - b := 0.5 * (a + !b / p.[1]) + b <- 0.5 * (a + b / p.[1]) quad () - k <- 4 + k <- 4 k <- k + 1 - if not !quadExecuted && p.[2] < 0.0 + if not quadExecuted && p.[2] < 0.0 then - b := sqrt !c - d <- !b + !b - a + b <- sqrt c + d <- b + b - a p.[1] <- 0.0 if d > 0.0 then p.[1] <- sqrt d - elif not !quadExecuted + elif not quadExecuted then if p.[1] > 0.0 then - b := (sqrt p.[2]) * 2.0 + p.[1] + b <- (sqrt p.[2]) * 2.0 + p.[1] else - b := -(sqrt p.[2]) * 2.0 + p.[1] + b <- -(sqrt p.[2]) * 2.0 + p.[1] - if !b <> 0.0 + if b <> 0.0 then p.[1] <- 0.0 else for k in 1..4 do r.[1, k] <- -e r.[2, k] <- 0.0 - quadExecuted := true + quadExecuted <- true quad () +// Return a tuple (area, x intersections, y intersections) +let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option = + let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha + let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha -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 + if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS then - -1.0 + None else - let phi_1 = phi_1 % Math.PI - let phi_2 = phi_2 % Math.PI + let phi_1 = phi_1 % Math.PI //(if phi_1 > Math.PI / 2.0 then phi_1 - Math.PI else phi_1) % Math.PI + let phi_2 = phi_2 % Math.PI //(if phi_2 > Math.PI / 2.0 then phi_2 - Math.PI else phi_2) % Math.PI let h2_tr, k2_tr, phi_2r = let cosphi = cos phi_1 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 @@ -562,7 +563,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = #if DEBUG_LOG for i in 0..4 do printf "cy[%d]=%f\n" i cy.[i] -#endif +#endif let py = Array.zeroCreate 5 let r = Array2D.zeroCreate 3 5 @@ -587,7 +588,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = py.[0] <- 1.0 cubicroots py r 3 - + elif abs cy.[2] > EPS then for i in 0..1 do @@ -609,15 +610,16 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = 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 + let ychk = Array.init nroots (fun _ -> Double.MaxValue) + let mutable nychk = 0 + for i in 1 .. nroots do + if abs r.[2, i] < EPS + then + ychk.[nychk] <- r.[1, i] * b1 + nychk <- nychk + 1 #if DEBUG_LOG - printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1 + printf "ROOT is Real, i=%d --> %f (B1=%f)\n" i r.[1, i] b1 #endif - |] Array.sortInPlace ychk #if DEBUG_LOG @@ -626,7 +628,6 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = printf "\t j=%d, ychk=%f\n" j ychk.[j] #endif - let nychk = Array.length ychk let mutable nintpts = 0 let xint = Array.zeroCreate 4 @@ -652,7 +653,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = #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 +#endif if abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS then @@ -661,7 +662,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = printf "first if x1. acc nintps=%d\n" nintpts #endif if nintpts > 4 - then + then returnValue <- -1.0 else xint.[nintpts-1] <- x1 @@ -677,7 +678,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = printf "first if x2. nintps=%d, Dx=%f (eps2=%f) \n" nintpts (abs (x2 - x1)) EPS #endif if nintpts > 4 - then + then returnValue <- -1.0 else xint.[nintpts-1] <- x2 @@ -692,26 +693,36 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : float = printf "i=%d, multiple roots: %f <--------> %f. continue\n" i ychk.[i] ychk.[i-1] #endif i <- i + 1 - - + + if returnValue = -1.0 then - returnValue + None else - 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 -> + let area = + 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 -> #if DEBUG_LOG - printf "one point is tangent\n" + 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 + 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 -> + | INTERSECTION_POINT -> #if DEBUG_LOG - printf "check twointpts\n" + 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 + 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 + if nintpts = 0 + then Some (float32 area, [||], [||]) + else + let xTransform : float32[] = Array.zeroCreate nintpts + let yTransform : float32[] = Array.zeroCreate nintpts + for i in 0 .. (nintpts - 1) do + xTransform.[i] <- float32 <| cos phi_1 * xint.[i] - sin phi_1 * yint.[i] + h1 + yTransform.[i] <- float32 <| sin phi_1 * xint.[i] + cos phi_1 * yint.[i] + k1 + Some (float32 area, xTransform, yTransform) \ No newline at end of file