From: Greg Burri Date: Wed, 9 Dec 2015 20:02:02 +0000 (+0100) Subject: First commit of the f# source code. X-Git-Tag: 1.0.11~86 X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=commitdiff_plain;h=ba64921fb9a0c36cd8cf802cbf1b2c0f79bc25f6 First commit of the f# source code. --- diff --git a/Parasitemia/Parasitemia.sln b/Parasitemia/Parasitemia.sln new file mode 100644 index 0000000..aeabd42 --- /dev/null +++ b/Parasitemia/Parasitemia.sln @@ -0,0 +1,28 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +# Visual Studio 14 +VisualStudioVersion = 14.0.22823.1 +MinimumVisualStudioVersion = 10.0.40219.1 +Project("{F2A71F9B-5D33-465A-A702-920D77279786}") = "Parasitemia", "Parasitemia\Parasitemia.fsproj", "{70838E65-F211-44FC-B28F-0ED1CA6E850F}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "WPF", "WPF\WPF.csproj", "{314FD78E-870E-4794-BB16-EA4586F2ABDB}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Any CPU = Debug|Any CPU + Release|Any CPU = Release|Any CPU + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {70838E65-F211-44FC-B28F-0ED1CA6E850F}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {70838E65-F211-44FC-B28F-0ED1CA6E850F}.Debug|Any CPU.Build.0 = Debug|Any CPU + {70838E65-F211-44FC-B28F-0ED1CA6E850F}.Release|Any CPU.ActiveCfg = Release|Any CPU + {70838E65-F211-44FC-B28F-0ED1CA6E850F}.Release|Any CPU.Build.0 = Release|Any CPU + {314FD78E-870E-4794-BB16-EA4586F2ABDB}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {314FD78E-870E-4794-BB16-EA4586F2ABDB}.Debug|Any CPU.Build.0 = Debug|Any CPU + {314FD78E-870E-4794-BB16-EA4586F2ABDB}.Release|Any CPU.ActiveCfg = Release|Any CPU + {314FD78E-870E-4794-BB16-EA4586F2ABDB}.Release|Any CPU.Build.0 = Release|Any CPU + EndGlobalSection + GlobalSection(SolutionProperties) = preSolution + HideSolutionNode = FALSE + EndGlobalSection +EndGlobal diff --git a/Parasitemia/Parasitemia/App.config b/Parasitemia/Parasitemia/App.config new file mode 100644 index 0000000..640d253 --- /dev/null +++ b/Parasitemia/Parasitemia/App.config @@ -0,0 +1,18 @@ + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/Parasitemia/Parasitemia/AssemblyInfo.fs b/Parasitemia/Parasitemia/AssemblyInfo.fs new file mode 100644 index 0000000..269e427 --- /dev/null +++ b/Parasitemia/Parasitemia/AssemblyInfo.fs @@ -0,0 +1,41 @@ +namespace Parasitemia.AssemblyInfo + +open System.Reflection +open System.Runtime.CompilerServices +open System.Runtime.InteropServices + +// General Information about an assembly is controlled through the following +// set of attributes. Change these attribute values to modify the information +// associated with an assembly. +[] +[] +[] +[] +[] +[] +[] +[] + +// Setting ComVisible to false makes the types in this assembly not visible +// to COM components. If you need to access a type in this assembly from +// COM, set the ComVisible attribute to true on that type. +[] + +// The following GUID is for the ID of the typelib if this project is exposed to COM +[] + +// Version information for an assembly consists of the following four values: +// +// Major Version +// Minor Version +// Build Number +// Revision +// +// You can specify all the values or you can default the Build and Revision Numbers +// by using the '*' as shown below: +// [] +[] +[] + +do + () \ No newline at end of file diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs new file mode 100644 index 0000000..b22a4fd --- /dev/null +++ b/Parasitemia/Parasitemia/Config.fs @@ -0,0 +1,9 @@ +module Config + +type Config = { + doGSigma1: float + doGSigma2: float + doGLowFreqPercentageReduction: float + + scale: float +} \ No newline at end of file diff --git a/Parasitemia/Parasitemia/EEOver.fs b/Parasitemia/Parasitemia/EEOver.fs new file mode 100644 index 0000000..e366b72 --- /dev/null +++ b/Parasitemia/Parasitemia/EEOver.fs @@ -0,0 +1,604 @@ +module EEOver + +open System + +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 a1b1 = a1 * b1 + let a2b2 = a2 * b2 + let area_1 = Math.PI * a1b1 + let area_2 = Math.PI * a2b2 + let relsize = a1b1 - a2b2 + + if relsize > 0.0 + then + if (h2_tr * h2_tr) / (a1 * a1) + (k2_tr * k2_tr) / (b1 * b1) < 1.0 + then area_2 + else 0.0 + + elif relsize < 0.0 + then + if ff < 0.0 + then area_1 + else 0.0 + + else + if abs (h1 - h2) < EPS && abs (k1 - k2) < EPS && abs (area_1 - area_2) < EPS + then area_1 + else 0.0 + +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 = + if abs x > a1 + then + if x < 0.0 then -a1 else a1 + else x + + let theta = + if y < 0.0 + then 2.0 * Math.PI - acos (x / a1) + else acos (x / a1) + + let eps_radian = 0.1 + + let x1 = a1 * cos (theta + eps_radian) + let y1 = b1 * sin (theta + eps_radian) + let x2 = a1 * cos (theta - eps_radian) + let y2 = b1 * sin (theta - eps_radian) + + let test1 = ellipse2tr x1 y1 aa bb cc dd ee ff + let test2 = ellipse2tr x2 y2 aa bb cc dd ee ff + + if test1 * test2 > 0.0 + then TANGENT_POINT + else INTERSECTION_POINT + + +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) = + if abs x.[0] > a1 + then x.[0] <- if x.[0] < 0.0 then -a1 else a1 + + let mutable theta1 = + if y.[0] < 0.0 + then 2.0 * Math.PI - acos (x.[0] / a1) + else acos (x.[0] / a1) + + if abs x.[1] > a1 + then x.[1] <- if x.[1] < 0.0 then -a1 else a1 + + let mutable theta2 = + if y.[1] < 0.0 + then 2.0 * Math.PI - acos (x.[1] / a1) + else acos (x.[1] / a1) + + if theta1 > theta2 + then + let tmp = theta1 + theta1 <- theta2 + theta2 <- tmp + + let xmid = a1 * cos ((theta1 + theta2) / 2.0) + let ymid = b1 * sin ((theta1 + theta2) / 2.0) + + if ellipse2tr xmid ymid aa bb cc dd ee ff > 0.0 + then + let tmp = theta1 + theta1 <- theta2 + theta2 <- tmp + + 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])) + + if area1 < 0.0 + then + area1 <- area1 + a1 * b1 + + let cosphi = cos (phi_1 - phi_2) + let sinphi = sin (phi_1 - phi_2) + + let mutable x1_tr = (x.[0] - h2_tr) * cosphi + (y.[0] - k2_tr) * -sinphi + let mutable y1_tr = (x.[0] - h2_tr) * sinphi + (y.[0] - k2_tr) * cosphi + let mutable x2_tr = (x.[1] - h2_tr) * cosphi + (y.[1] - k2_tr) * -sinphi + let mutable y2_tr = (x.[1] - h2_tr) * sinphi + (y.[1] - k2_tr) * cosphi + + if abs x1_tr > a2 + then + x1_tr <- if x1_tr < 0.0 then -a2 else a2 + + if y1_tr < 0.0 + then + theta1 <- 2.0 * Math.PI - acos (x1_tr / a2) + else + theta1 <- acos (x1_tr / a2) + + if abs x2_tr > a2 + then + x2_tr <- if x2_tr < 0.0 then -a2 else a2 + + if y2_tr < 0.0 + then + theta2 <- 2.0 * Math.PI - acos (x2_tr / a2) + else + theta2 <- acos (x2_tr / a2) + + if theta1 > theta2 + then + let tmp = theta1 + theta1 <- theta2 + theta2 <- tmp + + let xmid = a2 * cos ((theta1 + theta2) / 2.0) + let ymid = b2 * sin ((theta1 + theta2) / 2.0) + + let cosphi = cos (phi_2 - phi_1) + let sinphi = sin (phi_2 - phi_1) + let xmid_rt = xmid * cosphi + ymid * -sinphi + h2_tr + let ymid_rt = xmid * sinphi + ymid * cosphi + k2_tr + + if (xmid_rt * xmid_rt) / (a1 * a1) + (ymid_rt * ymid_rt) / (b1 * b1) > 1.0 + then + let tmp = theta1 + theta1 <- theta2 + theta2 <- tmp + + 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 area2 = 0.5 * (a2 * b2 * (theta2 - theta1) + trsign * abs (x1_tr * y2_tr - x2_tr * y1_tr)) + if area2 < 0.0 + then + area2 <- area2 + a2 * b2 + + 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 mutable tanpts = 0 + let mutable tanindex = 0 + 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; + + if tanpts <> 1 + then + -1.0 + else + match tanindex with + | 0 -> + xint.[0] <- xint.[2] + yint.[0] <- yint.[2] + | 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 + + +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 + let area_1 = Math.PI * a1b1 + let area_2 = Math.PI * a2b2 + + let theta = Array.zeroCreate 4 + + for i in 0..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) + + for j in 1..3 do + let tmp0 = theta.[j] + let tmp1 = xint.[j] + let tmp2 = yint.[j] + + let mutable k = j - 1 + while k >= 0 do + if theta.[k] <= tmp0 + then + k <- -1 + else + theta.[k+1] <- theta.[k] + xint.[k+1] <- xint.[k] + yint.[k+1] <- yint.[k] + k <- k - 1 + + theta.[k+1] <- tmp0 + xint.[k+1] <- tmp1 + yint.[k+1] <- tmp2 + + 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) + let sinphi = sin (phi_1 - phi_2) + + let theta_tr = Array.zeroCreate 4 + let xint_tr = Array.zeroCreate 4 + let yint_tr = Array.zeroCreate 4 + + for i in 0..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 + + 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) + let ymid = b1 * sin ((theta.[0] + theta.[1]) / 2.0) + + let mutable area2, area3, area4, area5 = 0.0, 0.0, 0.0, 0.0 + + 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])) + 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] + then + 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])) + else + area5 <- 0.5 * (a2b2 * (theta_tr.[0] - theta_tr.[3]) - abs (xint_tr.[3] * yint_tr.[0] - xint_tr.[0] * yint_tr.[3])) + else + area2 <- 0.5 * (a1b1 * (theta.[2] - theta.[1]) - abs (xint.[1] * yint.[2] - xint.[2] * yint.[1])) + 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 + area5 <- area5 + area_2 + + if area4 < 0.0 + then + area4 <- area4 + area_2 + + if area3 < 0.0 + then + area3 <- area3 + area_1 + + if area2 < 0.0 + then + area2 <- area2 + area_1 + + area1 + area2 + area3 + area4 + area5 + + +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 + + if d >= 0.0 + then + if b > 0.0 + then + b <- sqrt d + b + r.[1, 2] <- b + else + b <- -sqrt d + b + r.[1, 2] <- b + r.[1, 1] <- c / b + r.[2, 1] <- 0.0 + r.[2, 2] <- 0.0 + else + d <- sqrt -d + r.[2, 1] <- d + r.[2, 2] <- -d + r.[1, 1] <- b + r.[1, 2] <- b + +let private cubicroots (p: float[]) (r: float[,]) = + if p.[0] <> 1.0 then + for k in 1..3 do + p.[k] <- p.[k] / p.[0] + p.[0] <- 1.0 + let s = p.[1] / 3.0 + let mutable t = s * p.[1] + let mutable b = 0.5 * (s * (t / 1.5 - p.[2]) + p.[3]) + t <- (t - p.[2]) / 3.0 + let mutable c = t * t * t + let mutable d = b * b - c + + if d >= 0.0 + then + d <- ((sqrt d) + (abs b)) ** (1.0 / 3.0) + if d <> 0.0 + then + if b > 0.0 + then b <- -d + else b <- d + c <- t / b + d <- sqrt(0.75) * (b - c) + r.[2, 2] <- d + b <- b + c + c <- -0.5 * b - s + r.[1, 2] <- c + if b > 0.0 && s <= 0.0 || b < 0.0 && s > 0.0 + then + r.[1, 1] <- c + r.[2, 1] <- -d + r.[1, 3] <- b - s + r.[2, 3] <- 0.0 + else + r.[1, 1] <- b - s + r.[2, 1] <- 0.0 + r.[1, 3] <- c + r.[2, 3] <- -d + else + if b = 0.0 + then d <- (atan 1.0) / 1.5 + else d <- atan ((sqrt -d) / (abs b)) / 3.0 + + if b < 0.0 + then b <- 2.0 * (sqrt t) + else b <- -2.0 * (sqrt t) + + c <- (cos d) * b + t <- -(sqrt 0.75) * (sin d) * b - 0.5 * c + d <- -t - c - s + c <- c - s + t <- t - s + + if abs c > abs t + then + r.[1, 3] <- c + else + r.[1, 3] <- t + t <- c + + if abs d > abs t + then + r.[1, 2] <- d + else + r.[1, 2] <- t + t <- d + + 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 + then + for k in 1..4 do + 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 a = p.[2] - d + c := p.[4] + e * (e * a - p.[3]) + a <- a - d + + let quadExecuted = ref false + let quad () = + if not !quadExecuted + then + 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 + quadroots p r + for k in 1..4 do + r.[1,k] <- r.[1,k] - e + quadExecuted := true + + p.[1] <- 0.5 * a + p.[2] <- (p.[1] * p.[1] - !c) * 0.25 + p.[3] <- !b * !b / -64.0 + if p.[3] < 0.0 + then + cubicroots p r + let mutable k = 1 + while k < 4 do + if r.[2, k] = 0.0 && r.[1, k] > 0.0 + then + d <- r.[1, k] * 4.0 + a <- a + d + if a >= 0.0 && !b >= 0.0 + then + p.[1] <- sqrt d + elif a <= 0.0 && !b <= 0.0 + then + p.[1] <- sqrt d + else + p.[1] <- -(sqrt d) + b := 0.5 * (a + !b / p.[1]) + quad () + k <- 4 + k <- k + 1 + + if not !quadExecuted && p.[2] < 0.0 + then + b := sqrt !c + d <- !b + !b - a + p.[1] <- 0.0 + if d > 0.0 + then + p.[1] <- sqrt d + elif not !quadExecuted + then + if p.[1] > 0.0 + then + b := (sqrt p.[2]) * 2.0 + p.[1] + else + b := -(sqrt p.[2]) * 2.0 + p.[1] + + 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 + + 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 + + if a1 <= EPS || b1 <= EPS || a2 <= EPS || b2 <= EPS + then + -1.0 + else + 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 + (h2 - h1) * cosphi + (k2 - k1) * sinphi, (h1 - h2) * sinphi + (k2 - k1) * cosphi, (phi_2 - phi_1) % (2.0 * Math.PI) + + let cosphi = cos phi_2r + let cosphi2 = cosphi ** 2.0 + let sinphi = sin phi_2r + let sinphi2 = sinphi ** 2.0 + let cosphisinphi = 2.0 * cosphi * sinphi + let a22 = a2 ** 2.0 + let b22 = b2 ** 2.0 + let tmp0 = (cosphi * h2_tr + sinphi * k2_tr) / a22 + let tmp1 = (sinphi * h2_tr - cosphi * k2_tr) / b22 + let tmp2 = cosphi * h2_tr + sinphi * k2_tr + let tmp3 = sinphi * h2_tr - cosphi * k2_tr + + let aa = cosphi2 / a22 + sinphi2 / b22 + let bb = cosphisinphi / a22 - cosphisinphi / b22 + let cc = sinphi2 / a22 + cosphi2 / b22 + let dd = -2.0 * cosphi * tmp0 - 2.0 * sinphi * tmp1 + let ee = -2.0 * sinphi * tmp0 + 2.0 * cosphi * tmp1 + let ff = tmp2 * tmp2 / a22 + tmp3 * tmp3 / b22 - 1.0 + + let cy = [| + (a1 * (a1 * aa - dd) + ff) * (a1 * (a1 * aa + dd) + ff) + 2.0 * b1 * (a1 * a1 * (aa * ee - bb * dd) + ee * ff) + 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) + 2.0 * b1 * (b1 * b1 * cc * ee + a1 * a1 * (bb * dd - aa * ee)) + a1 * a1 * a1 * a1 * aa * aa + b1 * b1 * (a1 * a1 * (bb * bb - 2.0 * aa * cc) + b1 * b1 * cc * cc) + |] + + let py = Array.zeroCreate 5 + let r = Array2D.zeroCreate 3 5 + + let nroots = + if abs cy.[4] > EPS + then + for i in 0 .. 3 do + py.[4-i] <- cy.[i] / cy.[4] + py.[0] <- 1.0 + biquadroots py r + 4 + + elif abs cy.[3] > EPS + then + for i in 0..2 do + py.[3 - i] <- cy.[i] / cy.[3] + py.[0] <- 1.0 + cubicroots py r + 3 + + elif abs cy.[2] > EPS + then + for i in 0..1 do + py.[2-i] <- cy.[i] / cy.[2] + py.[0] <- 1.0 + quadroots py r + 2 + + elif abs cy.[1] > EPS + then + r.[1, 1] <- -cy.[0] / cy.[1] + r.[2, 1] <- 0.0 + 1 + + else + 0 + + let ychk = [| + for i in 1 .. nroots do + if abs r.[2, i] < EPS + then + yield r.[1, i] * b1 + |] + Array.sortInPlace ychk + + let nychk = Array.length ychk + let mutable nintpts = 0 + + let xint = Array.zeroCreate 4 + let yint = Array.zeroCreate 4 + + let mutable returnValue = 0.0 + + let mutable i = 0 + while returnValue = 0.0 && i < nychk do + if not (i < nychk - 1 && abs (ychk.[i] - ychk.[i+1]) < EPS / 2.0) + then + 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 abs (ellipse2tr x1 ychk.[i] aa bb cc dd ee ff) < EPS + then + nintpts <- nintpts + 1 + if nintpts > 4 + then + returnValue <- -1.0 + else + xint.[nintpts-1] <- x1 + yint.[nintpts-1] <- ychk.[i] + + 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 nintpts > 4 + then + returnValue <- -1.0 + else + xint.[nintpts-1] <- x2 + yint.[nintpts-1] <- ychk.[i] + i <- i + 1 + + + if returnValue = -1.0 + then + returnValue + 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 -> 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 + | 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 + diff --git a/Parasitemia/Parasitemia/Ellipse.fs b/Parasitemia/Parasitemia/Ellipse.fs new file mode 100644 index 0000000..62273b5 --- /dev/null +++ b/Parasitemia/Parasitemia/Ellipse.fs @@ -0,0 +1,269 @@ +module Ellipse + +open System +open System.Collections.Generic + +open Emgu.CV +open Emgu.CV.Structure + +open Utils +open MatchingEllipses + + +type private SearchExtremum = Minimum | Maximum + +let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) (xmax: float) (searchExtremum: SearchExtremum) : (float * float) = + let gr = 1.0 / 1.6180339887498948482 + let mutable a = xmin + let mutable b = xmax + let mutable c = b - gr * (b - a) + let mutable d = a + gr * (b - a) + + for i in 1 .. nbIter do + let mutable fc = f c + let mutable fd = f d + + if searchExtremum = Maximum + then + let tmp = fc + fc <- fd + fd <- tmp + + if fc < fd + then + b <- d; + d <- c; + c <- b - gr * (b - a); + else + a <- c; + c <- d; + d <- a + gr * (b - a); + + let x = (b + a) / 2.0 + x, f x + +let ellipse (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option = + let accuracy_extremum_search_1 = 4; + let accuracy_extremum_search_2 = 3; + + // p3 as the referencial. + let p1x = p1x - p3x + let p1y = p1y - p3y + + let p2x = p2x - p3x + let p2y = p2y - p3y + + // Convert to polar coordinates. + let alpha1 = atan m1 + let alpha2 = atan m2 + + let r1 = sqrt (p1x**2.0 + p1y**2.0) + let theta1 = atan2 p1y p1x + + let r2 = sqrt (p2x**2.0 + p2y**2.0) + let theta2 = atan2 p2y p2x + + let valid = + 4.0 * sin (alpha1 - theta1) * (-r1 * sin (alpha1 - theta1) + r2 * sin (alpha1 - theta2)) * + sin (alpha2 - theta2) * (-r1 * sin (alpha2 - theta1) + r2 * sin (alpha2 - theta2)) + + r1 * r2 * sin (alpha1 - alpha2) **2.0 * sin (theta1 - theta2) **2.0 < 0.0 + + if valid + then + let r theta = + (r1 * r2 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) / + (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0) + + let rabs = r >> abs + + // We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval. + let thetaTan, _ = goldenSectionSearch rabs accuracy_extremum_search_1 0.0 Math.PI Maximum + let rTan = r thetaTan + + let PTanx = rTan * cos thetaTan + let PTany = rTan * sin thetaTan + + let d1a = tan alpha1 + let d1b = -d1a * p1x + p1y + + let d2a = tan alpha2 + let d2b = -d2a * p2x + p2y + + let d3a = -1.0 / tan thetaTan + let d3b = -d3a * PTanx + PTany + + let Ux = -(d1b - d2b) / (d1a - d2a) + let Uy = -(d2a * d1b - d1a * d2b) / (d1a - d2a) + + let Vx = -(d1b - d3b) / (d1a - d3a) + let Vy = -(d3a * d1b - d1a * d3b) / (d1a - d3a) + + let Wx = p1x + (p2x - p1x) / 2.0 + let Wy = p1y + (p2y - p1y) / 2.0 + + let Zx = p1x + (PTanx - p1x) / 2.0 + let Zy = p1y + (PTany - p1y) / 2.0 + + let va = -(-Vy + Zy) / (Vx - Zx) + let vb = -(Zx * Vy - Vx * Zy) / (Vx - Zx) + + let ua = -(-Uy + Wy) / (Ux - Wx) + let ub = -(Wx * Uy - Ux * Wy) / (Ux - Wx) + + let cx = -(vb - ub) / (va - ua) + let cy = -(ua * vb - va * ub) / (va - ua) + + let rc = sqrt (cx**2.0 + cy**2.0) + let psi = atan2 cy cx + + let rellipse theta = + sqrt ( + rc**2.0 + (r1**2.0 * r2**2.0 * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2))**2.0 * sin (theta1 - theta2)**2.0) / + (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0)**2.0 - + (2.0 * r1 * r2 * rc * cos (theta - psi) * (r1 * (cos (alpha2 + theta - theta1 - theta2) - cos (alpha2 - theta) * cos (theta1 - theta2)) * sin (alpha1 - theta1) + r2 * (-cos (alpha1 + theta - theta1 - theta2) + cos (alpha1 - theta) * cos (theta1 - theta2)) * sin (alpha2 - theta2)) * sin (theta1 - theta2)) / + (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2))**2.0 - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2)**2.0)) + + // We search for an interval [theta_a, theta_b] and assume the function is unimodal in this interval. + let r1eTheta, r1e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.0 (Math.PI / 2.0) Maximum // Pi/2 and not pi because the period is Pi. + let r2eTheta, r2e = goldenSectionSearch rellipse accuracy_extremum_search_2 0.0 (Math.PI / 2.0) Minimum + + let rr1e = r r1eTheta + let r1ex = rr1e * cos r1eTheta + let r1ey = rr1e * sin r1eTheta + let mutable alpha = atan ((r1ey - cy) / (r1ex - cx)) + if alpha < 0.0 + then + alpha <- alpha + Math.PI + + // Ride off the p3 referential. + let cx = cx + p3x + let cy = cy + p3y + + Some { cx = cx; cy = cy; a = r1e; b = r2e; alpha = alpha } + else + None + + +let private vectorRotation (p1x: float) (p1y: float) (v1x: float) (v1y: float) (px: float) (py: float) : float = + let mutable rotation = 1.0 + if p1y > py + then + if v1x > 0.0 + then + rotation <- -1.0 + elif p1y < py + then + if v1x < 0.0 + then + rotation <- -1.0 + elif p1x > px + then + if v1y < 0.0 + then + rotation <- -1.0 + elif p1x < px + then + if v1y > 0.0 + then + rotation <- -1.0 + rotation + + +let private areVectorsValid (p1x: float) (p1y: float) (p2x: float) (p2y: float) (v1x: float) (v1y: float) (v2x: float) (v2y: float) : (float * float) option = + let m1 = -v1x / v1y + let m2 = -v2x / v2y + + let b1 = -m1 * p1x + p1y + let b2 = -m2 * p2x + p2y + let px = -((b1 - b2)/(m1 - m2)) + let py = -((m2 * b1 - m1 * b2)/(m1 - m2)) + + let rot1 = vectorRotation p1x p1y v1x v1y px py + let rot2 = vectorRotation p2x p2y v2x v2y px py + + if rot1 = rot2 || rot1 * atan2 (p1y - py) (p1x - px) + rot2 * atan2 (p2y - py) (p2x - px) <= 0.0 + then + None + else + Some (m1, m2) + + +let find (edges: Matrix) + (xDir: Image) + (yDir: Image) + (radiusRange: float * float) + (windowSize: int) + (factorNbPick: float) : Types.Ellipse list = + + let increment = windowSize / 4; + + let r1, r2 = radiusRange + let radiusTolerance = (r2 - r1) * 0.2 + + let minimumDistance = (r2 / 1.5) ** 2.0; + let squaredDistance x1 y1 x2 y2 = (x1 - x2) ** 2.0 + (y1 - y2) ** 2.0; + + let h = edges.Height + let w = edges.Width + + let mutable last_i, last_j = Int32.MaxValue, Int32.MaxValue + + let currentElements = List<(int * int)>() + + let edgesData = edges.Data + let xDirData = xDir.Data + let yDirData = yDir.Data + + let rng = Random() + + let ellipses = MatchingEllipses () + + for window_i in -windowSize + increment .. increment .. h - increment do + for window_j in -windowSize + increment .. increment .. w - increment do + + let window_i_begin = if window_i < 0 then 0 else window_i + let window_i_end = if window_i + windowSize - 1 >= h then h - 1 else window_i + windowSize - 1 + let window_j_begin = if window_j < 0 then 0 else window_j + let window_j_end = if window_j + windowSize - 1 >= w then w - 1 else window_j + windowSize - 1 + + // Remove old elements. + let indexFirstElement = currentElements.FindIndex(fun (_, pj) -> pj >= window_j) + if indexFirstElement > 0 + then currentElements.RemoveRange(0, indexFirstElement) + + // Add the new elements. + for j in window_j + windowSize - increment .. window_j + windowSize - 1 do + for i in window_i_begin .. window_i_end do + if j >= 0 && j < w && edgesData.[i, j] = 1uy + then currentElements.Add((i, j)) + + if currentElements.Count >= 10 + then + let mutable nbOfPicks = (float currentElements.Count) * factorNbPick |> int + while nbOfPicks > 0 do + let (p1y, p1x) as p1 = currentElements.[rng.Next(currentElements.Count)] + let (p2y, p2x) as p2 = currentElements.[rng.Next(currentElements.Count)] + let (p3y, p3x) as p3 = currentElements.[rng.Next(currentElements.Count)] + if p1 <> p2 && p1 <> p3 && p2 <> p3 + then + nbOfPicks <- nbOfPicks - 1 + let p1yf, p1xf = float p1y, float p1x + let p2yf, p2xf = float p2y, float p2x + let p3yf, p3xf = float p3y, float p3x + if squaredDistance p1xf p1yf p2xf p2yf >= minimumDistance && + squaredDistance p1xf p1yf p3xf p3yf >= minimumDistance && + squaredDistance p2xf p2yf p3xf p3yf >= minimumDistance + then + match areVectorsValid p1xf p1yf p2xf p2yf -xDirData.[p1y, p1x, 0] -yDirData.[p1y, p1x, 0] -xDirData.[p2y, p2x, 0] -yDirData.[p2y, p2x, 0] with + | Some (m1, m2) -> + match ellipse p1xf p1yf m1 p2xf p2yf m2 p3xf p3yf with + | Some e when e.cx > 0.0 && e.cx < (float w) - 1.0 && e.cy > 0.0 && e.cy < (float h) - 1.0 && + e.a >= r1 - radiusTolerance && e.a <= r2 + radiusTolerance && e.b >= r1 - radiusTolerance && e.b <= r2 + radiusTolerance -> + ellipses.Add e + | _ -> () + | _ -> () + + currentElements.Clear() + + ellipses.Ellipses + diff --git a/Parasitemia/Parasitemia/ImageAnalysis.fs b/Parasitemia/Parasitemia/ImageAnalysis.fs new file mode 100644 index 0000000..1215965 --- /dev/null +++ b/Parasitemia/Parasitemia/ImageAnalysis.fs @@ -0,0 +1,121 @@ +module ImageAnalysis + +open System +open System.Drawing + +open Emgu.CV +open Emgu.CV.Structure + +open Utils +open ImgTools +open Config +open Types + +type Result = { + RBCPositions : Point list + infectedRBCPositions : Point list + img: Image +} + +let doAnalysis (img: Image) (config: Config) : Result = + + let imgFloat = img.Convert() + use scaledImg = if config.scale = 1.0 then imgFloat else imgFloat.Resize(config.scale, CvEnum.Inter.Area) + + (*use scaledImg = + if config.scale = 1.0 + then + img + else + let m = new Mat() + CvInvoke.Resize(img, m, Size(roundInt (float img.Size.Width * config.scale), roundInt (float img.Size.Height * config.scale))) + m*) + + use green = scaledImg.Item(1) + + //use green = new Matrix(scaledImg.Size) + //CvInvoke.MixChannels(scaledImg, green, [| 1; 0 |]) + + //let greenMatrix = new Matrix(green.Height, green.Width, green.DataPointer) + + //let test = greenMatrix.[10, 10] + + use filteredGreen = (gaussianFilter green config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter green config.doGSigma2) + + use sobelKernel = + new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ] + [ 2.0f; 0.0f; -2.0f ] + [ 1.0f; 0.0f; -1.0f ]], Point(0, 0)) + + use xEdges = filteredGreen.Convolution(sobelKernel).Convert() + use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert() + + let xEdgesData = xEdges.Data + let yEdgesData = yEdges.Data + for r in 0..xEdges.Rows-1 do + xEdgesData.[r, 0, 0] <- 0.0 + xEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 + yEdgesData.[r, 0, 0] <- 0.0 + yEdgesData.[r, xEdges.Cols-1, 0] <- 0.0 + + for c in 0..xEdges.Cols-1 do + xEdgesData.[0, c, 0] <- 0.0 + xEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 + yEdgesData.[0, c, 0] <- 0.0 + yEdgesData.[xEdges.Rows-1, c, 0] <- 0.0 + + use magnitudes = new Matrix(xEdges.Size) + CvInvoke.CartToPolar(xEdges, yEdges, magnitudes, new Mat()) // Compute the magnitudes (without angles). + + let min = ref 0.0 + let minLocation = ref <| Point() + let max = ref 0.0 + let maxLocation = ref <| Point() + magnitudes.MinMax(min, max, minLocation, maxLocation) + + use magnitudesByte = ((magnitudes / !max) * 255.0).Convert() // Otsu from OpenCV only support 'byte'. + use edges = new Matrix(xEdges.Size) + let threshold = CvInvoke.Threshold(magnitudesByte, edges, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary) + thin edges + removeArea edges 12 + + saveMat (edges * 255.0) "edges.png" + + let radiusRange = config.scale * 20.0, config.scale * 40.0 + let windowSize = roundInt (1.6 * (snd radiusRange)) + let factorNbPick = 1.0; + let ellipses = Ellipse.find edges xEdges yEdges radiusRange windowSize factorNbPick + + drawEllipse img (List.head ellipses) (Bgr(0.0, 255.0, 255.0)) + saveImg img "ellipses.png" + + { RBCPositions = []; infectedRBCPositions = []; img = img } + + // + + (*use imageHSV = scaledImage.Convert() + let H, S = match imageHSV.Split() with // Warning: H is from 0 to 179°. + | [| H; S; _|] -> H, S + | _ -> failwith "unable to split the HSV channels" + + let hueShiftValue = 175 + // Modulo operator doesn't exist on matrix thus we have to apply a function to every pixels. + let correctedH : Image = H.Convert(fun b _ _ -> + (255 - int(b) * 255 / 179 + hueShiftValue) % 256 |> byte + ) + + let correctedS : Image = S.Not() + + let filteredH = correctedH.SmoothMedian(5) + let filteredS = correctedS.SmoothMedian(5)*) + + //let greenChannel = scaledImage.Item(1) + + //let filteredImage = (gaussianFilter greenChannel config.doGSigma1) - config.doGLowFreqPercentageReduction * (gaussianFilter greenChannel config.doGSigma2) + + // let filteredImage = greenChannel.ThresholdAdaptive(Gray(255.), CvEnum.AdaptiveThresholdType.GaussianC, CvEnum.ThresholdType.Binary, 61, Gray(5.0)) + // let thresholdedImage = filteredImage.CopyBlank() + + // CvInvoke.Threshold(filteredImage, thresholdedImage, 0., 255., CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.BinaryInv) |> ignore + + // filteredImage <| \ No newline at end of file diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs new file mode 100644 index 0000000..cc09405 --- /dev/null +++ b/Parasitemia/Parasitemia/ImgTools.fs @@ -0,0 +1,177 @@ +module ImgTools + +open System +open System.Drawing +open System.Collections.Generic + +open Emgu.CV +open Emgu.CV.Structure + +open Utils + +let gaussianFilter (img : Image<'TColor, 'TDepth>) (standardDeviation : float) : Image<'TColor, 'TDepth> = + let size = 2 * int (ceil (4.0 * standardDeviation)) + 1 + img.SmoothGaussian(size, size, standardDeviation, standardDeviation) + +// Zhang and Suen algorithm. +// Modify 'mat' in place. +let thin (mat: Matrix) = + let neighbors = [| + (-1, 0) // p2 + (-1, 1) // p3 + ( 0, 1) // p4 + ( 1, 1) // p5 + ( 1, 0) // p6 + ( 1, -1) // p7 + ( 0, -1) // p8 + (-1, -1) |] // p9 + + let w = mat.Width + let h = mat.Height + let mutable data1 = mat.Data + let mutable data2 = Array2D.zeroCreate h w + + // Return the list of neighbor values. + let neighborsValues (p1i, p1j) = + Array.map (fun (ni, nj) -> + let pi = p1i + ni + let pj = p1j + nj + if pi < 0 || pi >= h || pj < 0 || pj >= w then 0uy else data1.[pi, pj] + ) neighbors + + // Return the number of 01 pattern in 'values' in a circular way. + let pattern01 (values: byte[]) = + let mutable nb = 0 + let mutable lastValue = 255uy + for v in values do + if lastValue = 0uy && v = 1uy + then + nb <- nb + 1 + lastValue <- v + if lastValue = 0uy && values.[0] = 1uy + then + nb <- nb + 1 + nb + + let mutable pixelChanged = true + let mutable oddIteration = true + while pixelChanged do + pixelChanged <- false + for i in 0..h-1 do + for j in 0..w-1 do + if data1.[i, j] = 1uy + then + let values = neighborsValues (i, j) + let s = Array.reduce (+) values + if s >= 2uy && s <= 6uy && + pattern01 values = 1 && + (not oddIteration || (values.[0] * values.[2] * values.[4] = 0uy && values.[2] * values.[4] * values.[6] = 0uy)) && // Odd iteration. + (oddIteration || (values.[0] * values.[2] * values.[6] = 0uy && values.[0] * values.[4] * values.[6] = 0uy)) // Even iterations. + then + data2.[i, j] <- 0uy + pixelChanged <- true + else + data2.[i, j] <- 1uy + else + data2.[i, j] <- 0uy + + oddIteration <- not oddIteration + let tmp = data1 + data1 <- data2 + data2 <- tmp + + +// Remove all 8-connected pixels with an area equal or greater than 'areaSize'. +// Modify 'mat' in place. +let removeArea (mat: Matrix) (areaSize: int) = + let neighbors = [| + (-1, 0) // p2 + (-1, 1) // p3 + ( 0, 1) // p4 + ( 1, 1) // p5 + ( 1, 0) // p6 + ( 1, -1) // p7 + ( 0, -1) // p8 + (-1, -1) |] // p9 + + let mat' = new Matrix(mat.Size) + let w = mat'.Width + let h = mat'.Height + mat.CopyTo(mat') + + let data = mat.Data + let data' = mat'.Data + + for i in 0..h-1 do + for j in 0..w-1 do + if data'.[i, j] = 1uy + then + let neighborhood = List<(int*int)>() + let neighborsToCheck = List<(int*int)>() + neighborsToCheck.Add((i, j)) + data'.[i, j] <- 0uy + + let pop (l: List<'a>) : 'a = + let n = l.[l.Count - 1] + l.RemoveAt(l.Count - 1) + n + + while neighborsToCheck.Count > 0 do + let (ci, cj) = pop neighborsToCheck + neighborhood.Add((ci, cj)) + for (ni, nj) in neighbors do + let pi = ci + ni + let pj = cj + nj + if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy + then + neighborsToCheck.Add((pi, pj)) + data'.[pi, pj] <- 0uy + if neighborhood.Count <= areaSize + then + for (ni, nj) in neighborhood do + data.[ni, nj] <- 0uy + + +let saveImg (img: Image<'TColor, 'TDepth>) (name: string) = + img.Save("output/" + name) + + +let saveMat (mat: Matrix<'TDepth>) (name: string) = + use img = new Image(mat.Size) + mat.CopyTo(img) + saveImg img name + +(*let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = + let e' = Ellipse(PointF(float32 e.cx, float32 e.cy), SizeF(2.0f * float32 e.a, 2.0f * float32 e.b), float32 e.alpha) + img.Draw(e', color)*) + +let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) = + let x0, y0, x1, y1 = roundInt(x0), roundInt(y0), roundInt(x1), roundInt(y1) + + img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, 1); + +let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Types.Ellipse) (color: 'TColor) = + let cosAlpha = cos e.alpha + let sinAlpha = sin e.alpha + + let mutable x0 = 0.0 + let mutable y0 = 0.0 + let mutable first_iteration = true + + let n = 40 + let thetaIncrement = 2.0 * Math.PI / (float n) + + for theta in 0.0 .. thetaIncrement .. 2.0 * Math.PI do + let cosTheta = cos theta + let sinTheta = sin theta + let x = e.cx + cosAlpha * e.a * cosTheta - sinAlpha * e.b * sinTheta + let y = e.cy + sinAlpha * e.a * cosTheta + cosAlpha * e.b * sinTheta + + if not first_iteration + then + drawLine img color x0 y0 x y + else + first_iteration <- false + + x0 <- x + y0 <- y \ No newline at end of file diff --git a/Parasitemia/Parasitemia/KMedians.fs b/Parasitemia/Parasitemia/KMedians.fs new file mode 100644 index 0000000..5ba15db --- /dev/null +++ b/Parasitemia/Parasitemia/KMedians.fs @@ -0,0 +1,46 @@ +module KMedians + +open Emgu.CV +open Emgu.CV.Structure + +open System.Drawing + +let kmedians (mat: Matrix) (fgFactor: float) : Matrix = + let nbIteration = 3 + + let min = ref 0.0 + let minLocation = ref <| Point() + let max = ref 0.0 + let maxLocation = ref <| Point() + mat.MinMax(min, max, minLocation, maxLocation) + + let mutable bgValue = !max - (!max - !min) / 4.0 + let mutable fgValue = !min + (!max - !min) / 4.0 + let mutable d_bg = new Matrix(mat.Size) + let mutable d_fg = new Matrix(mat.Size) + + for i in 1..3 do + CvInvoke.Pow(mat - bgValue, 2.0, d_bg) + CvInvoke.Pow(mat - fgValue, 2.0, d_fg) + let fg = (d_fg * fgFactor).Cmp(d_bg, CvEnum.CmpType.LessThan) + + printfn "test" + + + (*d_bg = (imgFloat - color_bg) .^ 2.0; + d_fg = (imgFloat - color_fg) .^ 2.0; + fg = d_fg * Background_weight < d_bg; + imgFilteredFg = imgFloat; + imgFilteredFg(~fg) = nan; + color_fg = median(reshape(imgFilteredFg, [], 1), 'omitnan'); + imgFilteredBg = imgFloat; + imgFilteredBg(fg) = nan; + color_bg = median(reshape(imgFilteredBg, [], 1), 'omitnan'); + *) + + + new Matrix(mat.Size) + + + + diff --git a/Parasitemia/Parasitemia/MainWindow.xaml b/Parasitemia/Parasitemia/MainWindow.xaml new file mode 100644 index 0000000..4d1a8b0 --- /dev/null +++ b/Parasitemia/Parasitemia/MainWindow.xaml @@ -0,0 +1,30 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/Parasitemia/Parasitemia/MainWindow.xaml.fs b/Parasitemia/Parasitemia/MainWindow.xaml.fs new file mode 100644 index 0000000..26d30c9 --- /dev/null +++ b/Parasitemia/Parasitemia/MainWindow.xaml.fs @@ -0,0 +1,6 @@ +namespace Parasitemia.Views + +open FsXaml + +type MainWindow = XAML<"MainWindow.xaml"> + diff --git a/Parasitemia/Parasitemia/MatchingEllipses.fs b/Parasitemia/Parasitemia/MatchingEllipses.fs new file mode 100644 index 0000000..e4eae8d --- /dev/null +++ b/Parasitemia/Parasitemia/MatchingEllipses.fs @@ -0,0 +1,57 @@ +module MatchingEllipses + +open System +open System.Linq +open System.Collections.Generic + +open Types +open Utils + + +type EllipseScore = { mutable matchingScore: float; e: Ellipse } + +let matchingScoreThreshold1 = 0.5 +let matchingScoreThreshold2 = 2.0 + +let ellipseArea e = e.a * e.b * Math.PI + +type MatchingEllipses () = + let ellipses = List() + + member this.Add (e: Ellipse) = + + // dprintfn "new ellipse: %A, nb: %A" e ellipses.Count + + let areaE = ellipseArea e + + let mutable matchingScoreE = 0.0 + + for other in ellipses do + let areaOther = ellipseArea other.e + let commonArea = EEOver.EEOverlapArea e other.e + let matchingScore = (2.0 * commonArea / (areaE + areaOther)) ** 2.0 + if matchingScore >= matchingScoreThreshold1 + then + other.matchingScore <- other.matchingScore + matchingScore + matchingScoreE <- matchingScoreE + matchingScore + printfn "Score" + + ellipses.Add({ matchingScore = matchingScoreE; e = e }) + + member this.Ellipses : Ellipse list = + dprintfn "Number of ellipse: %A" ellipses.Count + + (*let sortedEllipses = + List.filter (fun e -> e.matchingScore >= matchingScoreThreshold2) (List.ofSeq ellipses) + |> List.sortWith (fun e1 e2 -> e2.matchingScore.CompareTo(e1.matchingScore))*) + + List.filter (fun e -> e.matchingScore >= matchingScoreThreshold2) (List.ofSeq ellipses) + |> List.sortWith (fun e1 e2 -> e2.matchingScore.CompareTo(e1.matchingScore)) + |> List.map (fun { e = e } -> e) + + // ellipses.Where(fun e -> e.matchingScore >= matchingScoreThreshold2) + // ([], fun acc { matchingScore = score; e = e } -> if score >= matchingScoreThreshold2 then e :: acc else acc) + + + + diff --git a/Parasitemia/Parasitemia/NumericUpDown.xaml b/Parasitemia/Parasitemia/NumericUpDown.xaml new file mode 100644 index 0000000..ae20a5b --- /dev/null +++ b/Parasitemia/Parasitemia/NumericUpDown.xaml @@ -0,0 +1,21 @@ + + + + + + + + + + +