X-Git-Url: http://git.euphorik.ch/?a=blobdiff_plain;f=Parasitemia%2FParasitemiaCore%2FEEOver.fs;h=9a77f48ab44d95f6fc5a0787d89c8abc06c6030f;hb=1b8e45987bde692ab5602c281f878707f70459b7;hp=a1358e99cf4a027e87e66dc0de9712d012685730;hpb=4bfa3cbdc6145e6944f02e24829ab2ef3a851ac1;p=master-thesis.git diff --git a/Parasitemia/ParasitemiaCore/EEOver.fs b/Parasitemia/ParasitemiaCore/EEOver.fs index a1358e9..9a77f48 100644 --- a/Parasitemia/ParasitemiaCore/EEOver.fs +++ b/Parasitemia/ParasitemiaCore/EEOver.fs @@ -1,8 +1,9 @@ -module ParasitemiaCore.EEOver +// Translation from https://github.com/chraibi/EEOver. +module ParasitemiaCore.EEOver open System -let private EPS = 1.0e-5 +let private EPS = 1.0e-7 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 @@ -183,7 +184,7 @@ let private twointpts (x: float[]) (y: float[]) (a1: float) (b1: float) (phi_1: 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 + for i = 0 to 2 do if istanpt xint.[i] yint.[i] a1 b2 aa bb cc dd ee ff = TANGENT_POINT then tanpts <- tanpts + 1 @@ -207,7 +208,6 @@ let private threeintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) () twointpts xint yint a1 b1 phi_1 a2 b2 h2_tr k2_tr phi_2 aa bb cc dd ee ff - 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 = let a1b1 = a1 * b1 let a2b2 = a2 * b2 @@ -216,18 +216,18 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( let theta = Array.zeroCreate 4 - for i in 0 .. 3 do + for i = 0 to 3 do if abs xint.[i] > a1 then 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 + for k = 0 to 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 + for j = 1 to 3 do let tmp0 = theta.[j] let tmp1 = xint.[j] let tmp2 = yint.[j] @@ -253,7 +253,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( #if DEBUG_LOG printf "AFTER sorting\n" - for k in 0..3 do + for k = 0 to 3 do printf "k=%d: Theta = %f, xint=%f, yint=%f\n" k theta.[k] xint.[k] yint.[k] #endif @@ -266,7 +266,7 @@ let private fourintpts (xint: float[]) (yint: float[]) (a1: float) (b1: float) ( let xint_tr = Array.zeroCreate 4 let yint_tr = Array.zeroCreate 4 - for i in 0..3 do + for i = 0 to 3 do xint_tr.[i] <- (xint.[i] - h2_tr) * cosphi + (yint.[i] - k2_tr) * -sinphi yint_tr.[i] <- (xint.[i] - h2_tr) * sinphi + (yint.[i] - k2_tr) * cosphi @@ -301,28 +301,28 @@ 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 + printf "\n\t\t-------------> area5 is negative (%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 + printf "\n\t\t-------------> area4 is negative (%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 + printf "\n\t\t-------------> area3 is negative (%f). Add: pi*A2*B2=%f <------------\n" area3 area_1 #endif area3 <- area3 + area_1 if area2 < 0.0 then #if DEBUG_LOG - printf "\n\t\t-------------> area2 is negativ (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1 + printf "\n\t\t-------------> area2 is negative (%f). Add: pi*A2*B2=%f <------------\n" area2 area_1 #endif area2 <- area2 + area_1 @@ -358,7 +358,7 @@ let private quadroots (p: float[]) (r: float[,]) = let private cubicroots (p: float[]) (r: float[,]) = if p.[0] <> 1.0 then - for k in 1..3 do + for k = 1 to 3 do p.[k] <- p.[k] / p.[0] p.[0] <- 1.0 let s = p.[1] / 3.0 @@ -423,13 +423,13 @@ let private cubicroots (p: float[]) (r: float[,]) = t <- d r.[1, 1] <- t - for k in 1..3 do + for k = 1 to 3 do r.[2, k] <- 0.0 let private biquadroots (p: float[]) (r: float[,]) = if p.[0] <> 1.0 then - for k in 1..4 do + for k = 1 to 4 do p.[k] <- p.[k] / p.[0] p.[0] <- 1.0 let e = 0.25 * p.[1] @@ -447,13 +447,13 @@ let private biquadroots (p: float[]) (r: float[,]) = then p.[2] <- c / b quadroots p r - for k in 1..2 do - for j in 1..2 do + for k = 1 to 2 do + for j = 1 to 2 do r.[j, k+2] <- r.[j, k] p.[1] <- -p.[1] p.[2] <- b quadroots p r - for k in 1..4 do + for k = 1 to 4 do r.[1,k] <- r.[1,k] - e quadExecuted <- true @@ -502,14 +502,16 @@ let private biquadroots (p: float[]) (r: float[,]) = then p.[1] <- 0.0 else - for k in 1..4 do + for k = 1 to 4 do r.[1, k] <- -e r.[2, k] <- 0.0 quadExecuted <- true quad () -// Return a tuple (area, x intersections, y intersections) +/// +/// 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 @@ -518,8 +520,8 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] then None else - 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 phi_1 = phi_1 % Math.PI + let phi_2 = phi_2 % Math.PI let h2_tr, k2_tr, phi_2r = let cosphi = cos phi_1 let sinphi = sin phi_1 @@ -557,7 +559,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] |] #if DEBUG_LOG - for i in 0..4 do + for i = 0 to 4 do printf "cy[%d]=%f\n" i cy.[i] #endif @@ -567,11 +569,11 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] let nroots = if abs cy.[4] > EPS then - for i in 0 .. 3 do + for i = 0 to 3 do py.[4-i] <- cy.[i] / cy.[4] py.[0] <- 1.0 #if DEBUG_LOG - for i in 0..4 do + for i to 0 to 4 do printf "py[%d]=%f\n" i py.[i] #endif biquadroots py r @@ -579,7 +581,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] elif abs cy.[3] > EPS then - for i in 0..2 do + for i = 0 to 2 do py.[3 - i] <- cy.[i] / cy.[3] py.[0] <- 1.0 cubicroots py r @@ -587,7 +589,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] elif abs cy.[2] > EPS then - for i in 0..1 do + for i = 0 to 1 do py.[2-i] <- cy.[i] / cy.[2] py.[0] <- 1.0 quadroots py r @@ -608,7 +610,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] let ychk = Array.init nroots (fun _ -> Double.MaxValue) let mutable nychk = 0 - for i in 1 .. nroots do + for i = 1 to nroots do if abs r.[2, i] < EPS then ychk.[nychk] <- r.[1, i] * b1 @@ -620,7 +622,7 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] #if DEBUG_LOG printf "nychk=%d\n" ychk.Length - for j in 0 .. ychk.Length - 1 do + for j = 0 to ychk.Length - 1 do printf "\t j=%d, ychk=%f\n" j ychk.[j] #endif @@ -647,8 +649,8 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] 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) + 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) + 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 @@ -690,7 +692,6 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] #endif i <- i + 1 - if returnValue = -1.0 then None @@ -713,8 +714,13 @@ let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] | 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, [||], [||]) + + if area = -1.0 + then + None + elif nintpts = 0 + then + Some (float32 area, [||], [||]) else let xTransform : float32[] = Array.zeroCreate nintpts let yTransform : float32[] = Array.zeroCreate nintpts