First commit of the f# source code.
authorGreg Burri <greg.burri@gmail.com>
Wed, 9 Dec 2015 20:02:02 +0000 (21:02 +0100)
committerGreg Burri <greg.burri@gmail.com>
Wed, 9 Dec 2015 20:02:02 +0000 (21:02 +0100)
22 files changed:
Parasitemia/Parasitemia.sln [new file with mode: 0644]
Parasitemia/Parasitemia/App.config [new file with mode: 0644]
Parasitemia/Parasitemia/AssemblyInfo.fs [new file with mode: 0644]
Parasitemia/Parasitemia/Config.fs [new file with mode: 0644]
Parasitemia/Parasitemia/EEOver.fs [new file with mode: 0644]
Parasitemia/Parasitemia/Ellipse.fs [new file with mode: 0644]
Parasitemia/Parasitemia/ImageAnalysis.fs [new file with mode: 0644]
Parasitemia/Parasitemia/ImgTools.fs [new file with mode: 0644]
Parasitemia/Parasitemia/KMedians.fs [new file with mode: 0644]
Parasitemia/Parasitemia/MainWindow.xaml [new file with mode: 0644]
Parasitemia/Parasitemia/MainWindow.xaml.fs [new file with mode: 0644]
Parasitemia/Parasitemia/MatchingEllipses.fs [new file with mode: 0644]
Parasitemia/Parasitemia/NumericUpDown.xaml [new file with mode: 0644]
Parasitemia/Parasitemia/NumericUpDown.xaml.fs [new file with mode: 0644]
Parasitemia/Parasitemia/Parasitemia.fsproj [new file with mode: 0644]
Parasitemia/Parasitemia/Program.fs [new file with mode: 0644]
Parasitemia/Parasitemia/Types.fs [new file with mode: 0644]
Parasitemia/Parasitemia/Utils.fs [new file with mode: 0644]
Parasitemia/Parasitemia/packages.config [new file with mode: 0644]
Parasitemia/WPF/BitmapSourceConverter.cs [new file with mode: 0644]
Parasitemia/WPF/Properties/AssemblyInfo.cs [new file with mode: 0644]
Parasitemia/WPF/WPF.csproj [new file with mode: 0644]

diff --git a/Parasitemia/Parasitemia.sln b/Parasitemia/Parasitemia.sln
new file mode 100644 (file)
index 0000000..aeabd42
--- /dev/null
@@ -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 (file)
index 0000000..640d253
--- /dev/null
@@ -0,0 +1,18 @@
+<?xml version="1.0" encoding="utf-8"?>
+<configuration>
+  <startup>
+    <supportedRuntime version="v4.0" sku=".NETFramework,Version=v4.6" />
+  </startup>
+  <runtime>
+    <assemblyBinding xmlns="urn:schemas-microsoft-com:asm.v1">
+      <dependentAssembly>
+        <assemblyIdentity name="FSharp.Core" publicKeyToken="b03f5f7f11d50a3a" culture="neutral" />
+        <bindingRedirect oldVersion="0.0.0.0-4.4.0.0" newVersion="4.4.0.0" />
+      </dependentAssembly>
+      <dependentAssembly>
+        <assemblyIdentity name="Castle.Core" publicKeyToken="407dd0808d44fbdc" culture="neutral" />
+        <bindingRedirect oldVersion="0.0.0.0-1.1.0.0" newVersion="1.1.0.0" />
+      </dependentAssembly>
+    </assemblyBinding>
+  </runtime>
+</configuration>
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/AssemblyInfo.fs b/Parasitemia/Parasitemia/AssemblyInfo.fs
new file mode 100644 (file)
index 0000000..269e427
--- /dev/null
@@ -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.
+[<assembly: AssemblyTitle("Parasitemia")>]
+[<assembly: AssemblyDescription("")>]
+[<assembly: AssemblyConfiguration("")>]
+[<assembly: AssemblyCompany("")>]
+[<assembly: AssemblyProduct("Parasitemia")>]
+[<assembly: AssemblyCopyright("Copyright ©  2015")>]
+[<assembly: AssemblyTrademark("")>]
+[<assembly: AssemblyCulture("")>]
+
+// 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.
+[<assembly: ComVisible(false)>]
+
+// The following GUID is for the ID of the typelib if this project is exposed to COM
+[<assembly: Guid("70838e65-f211-44fc-b28f-0ed1ca6e850f")>]
+
+// 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:
+// [<assembly: AssemblyVersion("1.0.*")>]
+[<assembly: AssemblyVersion("1.0.0.0")>]
+[<assembly: AssemblyFileVersion("1.0.0.0")>]
+
+do
+    ()
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/Config.fs b/Parasitemia/Parasitemia/Config.fs
new file mode 100644 (file)
index 0000000..b22a4fd
--- /dev/null
@@ -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 (file)
index 0000000..e366b72
--- /dev/null
@@ -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<float> 5
+        let r = Array2D.zeroCreate<float> 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 (file)
index 0000000..62273b5
--- /dev/null
@@ -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<byte>)
+         (xDir: Image<Gray, float>) 
+         (yDir: Image<Gray, float>) 
+         (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 (file)
index 0000000..1215965
--- /dev/null
@@ -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<Bgr, byte>
+}
+
+let doAnalysis (img: Image<Bgr, byte>) (config: Config) : Result =
+
+    let imgFloat = img.Convert<Bgr, float32>()
+    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<byte>(scaledImg.Size)
+    //CvInvoke.MixChannels(scaledImg, green, [| 1; 0 |])
+
+    //let greenMatrix = new Matrix<byte>(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<Gray, float>()
+    use yEdges = filteredGreen.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+            
+    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<float>(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<byte>() // Otsu from OpenCV only support 'byte'.
+    use edges = new Matrix<byte>(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<Hsv, uint8>()
+    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<Gray, byte> = H.Convert(fun b _ _ ->
+        (255 - int(b) * 255 / 179 + hueShiftValue) % 256 |> byte
+    )
+   
+    let correctedS : Image<Gray, byte> = 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 (file)
index 0000000..cc09405
--- /dev/null
@@ -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<byte>) =    
+    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<byte> 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<byte>) (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<byte>(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<Gray, 'TDeph>(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 (file)
index 0000000..5ba15db
--- /dev/null
@@ -0,0 +1,46 @@
+module KMedians
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+open System.Drawing
+
+let kmedians (mat: Matrix<float32>) (fgFactor: float) : Matrix<byte> =
+    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<float32>(mat.Size)
+    let mutable d_fg = new Matrix<float32>(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<byte>(mat.Size)
+
+
+
+
diff --git a/Parasitemia/Parasitemia/MainWindow.xaml b/Parasitemia/Parasitemia/MainWindow.xaml
new file mode 100644 (file)
index 0000000..4d1a8b0
--- /dev/null
@@ -0,0 +1,30 @@
+<Window xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
+        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:d="http://schemas.microsoft.com/expression/blend/2008" xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006" mc:Ignorable="d" x:Name="MainWindow" Height="931.638" Width="1360.61">
+   <Grid x:Name="mainGrid" Margin="10" RenderTransformOrigin="0.491,0.524">
+      <Grid.RowDefinitions>
+         <RowDefinition/>
+         <RowDefinition Height="100"/>
+      </Grid.RowDefinitions>
+      <Grid.ColumnDefinitions>
+         <ColumnDefinition/>
+         <ColumnDefinition Width="300"/>
+      </Grid.ColumnDefinitions>
+      <Image x:Name="img" HorizontalAlignment="Stretch" Margin="10" VerticalAlignment="Stretch"/>
+      <TextBlock x:Name="txtLog" Margin="10" TextWrapping="Wrap" FontFamily="Monaco" Grid.Row="2" Background="#FF95E633" Grid.ColumnSpan="2"/>
+      <Grid Grid.Column="1" Margin="10" >
+         <Grid.RowDefinitions>
+            <RowDefinition Height="40"/>
+            <RowDefinition Height="40"/>
+            <RowDefinition Height="40"/>
+            <RowDefinition />
+         </Grid.RowDefinitions>
+         <Grid.ColumnDefinitions>
+            <ColumnDefinition/>
+            <ColumnDefinition/>
+         </Grid.ColumnDefinitions>
+         <Label Content="Ecart type 1" Margin="0,7" VerticalAlignment="Center" HorizontalAlignment="Left" HorizontalContentAlignment="Stretch" />
+         <Label Content="Facteur deuxième gaussienne" Margin="0,7" VerticalAlignment="Center" HorizontalAlignment="Left" HorizontalContentAlignment="Stretch" Grid.Row="1" />
+         <Label Content="Facteur deuxième gaussienne" Margin="0,7" VerticalAlignment="Center" HorizontalAlignment="Left" HorizontalContentAlignment="Stretch" Grid.Row="2" />
+      </Grid>
+   </Grid>
+</Window>
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/MainWindow.xaml.fs b/Parasitemia/Parasitemia/MainWindow.xaml.fs
new file mode 100644 (file)
index 0000000..26d30c9
--- /dev/null
@@ -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 (file)
index 0000000..e4eae8d
--- /dev/null
@@ -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<EllipseScore>()
+        
+    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 (file)
index 0000000..ae20a5b
--- /dev/null
@@ -0,0 +1,21 @@
+<UserControl
+               xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
+               xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
+               xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
+      xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
+      mc:Ignorable="d" d:DesignWidth="259.5" d:DesignHeight="45.5"
+               >
+   <Grid>
+      <Grid.RowDefinitions>
+         <RowDefinition />
+         <RowDefinition />
+      </Grid.RowDefinitions>
+      <Grid.ColumnDefinitions>
+         <ColumnDefinition Width="*" />
+         <ColumnDefinition Width="Auto" />
+      </Grid.ColumnDefinitions>
+      <Button x:Name="upButton" Grid.Column="1" Content="^"/>
+      <Button x:Name="downButton" Grid.Column="1" Grid.Row="1" Content="v"/>
+      <TextBox x:Name="input" Grid.RowSpan="2" />
+   </Grid>
+</UserControl>
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/NumericUpDown.xaml.fs b/Parasitemia/Parasitemia/NumericUpDown.xaml.fs
new file mode 100644 (file)
index 0000000..b30e59a
--- /dev/null
@@ -0,0 +1,17 @@
+namespace Parasitemia.Views
+
+open System
+open System.Windows
+open System.Windows.Data
+open System.Windows.Input
+
+open FsXaml
+
+type NumericUpDown = XAML<"NumericUpDown.xaml", true>
+    
+type NumericUpDownEvents = Up | Down
+
+type NumericUpDownController() = 
+    inherit UserControlViewController<NumericUpDown>()
+
+
diff --git a/Parasitemia/Parasitemia/Parasitemia.fsproj b/Parasitemia/Parasitemia/Parasitemia.fsproj
new file mode 100644 (file)
index 0000000..d3a8fe4
--- /dev/null
@@ -0,0 +1,139 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="14.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <Import Project="$(MSBuildExtensionsPath)\$(MSBuildToolsVersion)\Microsoft.Common.props" Condition="Exists('$(MSBuildExtensionsPath)\$(MSBuildToolsVersion)\Microsoft.Common.props')" />
+  <PropertyGroup>
+    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
+    <Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
+    <SchemaVersion>2.0</SchemaVersion>
+    <ProjectGuid>70838e65-f211-44fc-b28f-0ed1ca6e850f</ProjectGuid>
+    <OutputType>WinExe</OutputType>
+    <RootNamespace>Parasitemia</RootNamespace>
+    <AssemblyName>Parasitemia</AssemblyName>
+    <TargetFrameworkVersion>v4.6</TargetFrameworkVersion>
+    <AutoGenerateBindingRedirects>true</AutoGenerateBindingRedirects>
+    <TargetFSharpCoreVersion>4.4.0.0</TargetFSharpCoreVersion>
+    <Name>Parasitemia</Name>
+    <NuGetPackageImportStamp>
+    </NuGetPackageImportStamp>
+    <TargetFrameworkProfile />
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
+    <DebugSymbols>true</DebugSymbols>
+    <DebugType>full</DebugType>
+    <Optimize>false</Optimize>
+    <Tailcalls>false</Tailcalls>
+    <OutputPath>bin\Debug\</OutputPath>
+    <DefineConstants>DEBUG;TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+    <PlatformTarget>x64</PlatformTarget>
+    <DocumentationFile>bin\Debug\Parasitemia.XML</DocumentationFile>
+    <Prefer32Bit>false</Prefer32Bit>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <Tailcalls>true</Tailcalls>
+    <OutputPath>bin\Release\</OutputPath>
+    <DefineConstants>TRACE</DefineConstants>
+    <WarningLevel>3</WarningLevel>
+    <PlatformTarget>x64</PlatformTarget>
+    <DocumentationFile>bin\Release\Parasitemia.XML</DocumentationFile>
+    <Prefer32Bit>false</Prefer32Bit>
+  </PropertyGroup>
+  <PropertyGroup>
+    <MinimumVisualStudioVersion Condition="'$(MinimumVisualStudioVersion)' == ''">11</MinimumVisualStudioVersion>
+  </PropertyGroup>
+  <Choose>
+    <When Condition="'$(VisualStudioVersion)' == '11.0'">
+      <PropertyGroup Condition="Exists('$(MSBuildExtensionsPath32)\..\Microsoft SDKs\F#\3.0\Framework\v4.0\Microsoft.FSharp.Targets')">
+        <FSharpTargetsPath>$(MSBuildExtensionsPath32)\..\Microsoft SDKs\F#\3.0\Framework\v4.0\Microsoft.FSharp.Targets</FSharpTargetsPath>
+      </PropertyGroup>
+    </When>
+    <Otherwise>
+      <PropertyGroup Condition="Exists('$(MSBuildExtensionsPath32)\Microsoft\VisualStudio\v$(VisualStudioVersion)\FSharp\Microsoft.FSharp.Targets')">
+        <FSharpTargetsPath>$(MSBuildExtensionsPath32)\Microsoft\VisualStudio\v$(VisualStudioVersion)\FSharp\Microsoft.FSharp.Targets</FSharpTargetsPath>
+      </PropertyGroup>
+    </Otherwise>
+  </Choose>
+  <Import Project="$(FSharpTargetsPath)" />
+  <ItemGroup>
+    <Compile Include="AssemblyInfo.fs" />
+    <None Include="NumericUpDown.xaml" />
+    <Compile Include="NumericUpDown.xaml.fs" />
+    <Resource Include="MainWindow.xaml" />
+    <Compile Include="MainWindow.xaml.fs" />
+    <Compile Include="Utils.fs" />
+    <Compile Include="Types.fs" />
+    <Compile Include="ImgTools.fs" />
+    <Compile Include="Config.fs" />
+    <Compile Include="EEOver.fs" />
+    <Compile Include="MatchingEllipses.fs" />
+    <Compile Include="Ellipse.fs" />
+    <Compile Include="KMedians.fs" />
+    <Compile Include="ImageAnalysis.fs" />
+    <Compile Include="Program.fs" />
+    <None Include="App.config" />
+    <Content Include="packages.config" />
+  </ItemGroup>
+  <ItemGroup>
+    <Reference Include="Emgu.CV">
+      <HintPath>..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.CV.dll</HintPath>
+    </Reference>
+    <Reference Include="Emgu.Util">
+      <HintPath>..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.Util.dll</HintPath>
+    </Reference>
+    <Reference Include="FSharp.Data.TypeProviders" />
+    <Reference Include="FSharp.ViewModule.Core.Wpf">
+      <HintPath>..\packages\FSharp.ViewModule.Core.0.9.9.1\lib\net45\FSharp.ViewModule.Core.Wpf.dll</HintPath>
+      <Private>True</Private>
+    </Reference>
+    <Reference Include="FsXaml.Wpf">
+      <HintPath>..\packages\FsXaml.Wpf.0.9.9\lib\net45\FsXaml.Wpf.dll</HintPath>
+      <Private>True</Private>
+    </Reference>
+    <Reference Include="FsXaml.Wpf.TypeProvider">
+      <HintPath>..\packages\FsXaml.Wpf.0.9.9\lib\net45\FsXaml.Wpf.TypeProvider.dll</HintPath>
+      <Private>True</Private>
+    </Reference>
+    <Reference Include="log4net">
+      <HintPath>..\packages\log4net.1.2.10\lib\2.0\log4net.dll</HintPath>
+      <Private>True</Private>
+    </Reference>
+    <Reference Include="mscorlib" />
+    <Reference Include="FSharp.Core, Version=$(TargetFSharpCoreVersion), Culture=neutral, PublicKeyToken=b03f5f7f11d50a3a">
+      <Private>True</Private>
+    </Reference>
+    <Reference Include="PresentationCore" />
+    <Reference Include="PresentationFramework" />
+    <Reference Include="System" />
+    <Reference Include="System.Core" />
+    <Reference Include="System.Data" />
+    <Reference Include="System.Data.Linq" />
+    <Reference Include="System.Drawing" />
+    <Reference Include="System.Numerics" />
+    <Reference Include="System.Windows.Interactivity">
+      <HintPath>..\packages\Expression.Blend.Sdk.1.0.2\lib\net45\System.Windows.Interactivity.dll</HintPath>
+      <Private>True</Private>
+    </Reference>
+    <Reference Include="System.Xaml" />
+    <Reference Include="System.Xml" />
+    <Reference Include="WindowsBase" />
+  </ItemGroup>
+  <ItemGroup>
+    <ProjectReference Include="..\WPF\WPF.csproj">
+      <Name>WPF</Name>
+      <Project>{314fd78e-870e-4794-bb16-ea4586f2abdb}</Project>
+      <Private>True</Private>
+    </ProjectReference>
+  </ItemGroup>
+  <PropertyGroup>
+    <PostBuildEvent>xcopy "D:\Emgu\emgucv-windows-universal 3.0.0.2157\bin\x64\*" "$(TargetDir)" /Y /D</PostBuildEvent>
+  </PropertyGroup>
+  <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
+       Other similar extension points exist, see Microsoft.Common.targets.
+  <Target Name="BeforeBuild">
+  </Target>
+  <Target Name="AfterBuild">
+  </Target>
+  -->
+</Project>
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/Program.fs b/Parasitemia/Parasitemia/Program.fs
new file mode 100644 (file)
index 0000000..56672bd
--- /dev/null
@@ -0,0 +1,52 @@
+module Parasitemia.Main
+
+open System.IO
+open System.Windows
+open System.Windows.Media
+open System.Windows.Markup
+open System.Windows.Shapes
+open System.Windows.Controls
+open System.Drawing
+open System.Diagnostics
+
+open Emgu.CV
+open Emgu.CV.Structure
+open Emgu.CV.WPF
+
+open Config
+
+let display (window : Views.MainWindow) (img : IImage) = 
+   let imgControl = window.Root.FindName("img") :?> Controls.Image
+   imgControl.Source <- BitmapSourceConvert.ToBitmapSource(img)
+
+let log (window : Views.MainWindow) (mess : string) =
+   let txtLog = window.Root.FindName("txtLog") :?> Controls.TextBlock
+   txtLog.Text <- txtLog.Text + mess + "\n"
+   
+[<System.STAThreadAttribute>]
+do
+    Utils.dprintfn "OpenCV test"
+
+    let app = new Application()
+    let mainWindow = Views.MainWindow()
+
+    Utils.log <- (fun m -> log mainWindow m)
+
+    let config = {
+        doGSigma1 = 1.5
+        doGSigma2 = 20.0
+        doGLowFreqPercentageReduction = 0.9
+        scale = 1.0
+    }
+   
+    //use img = new Image<Rgb, byte>("../../../../src/Tests_hough/images/1508133543-7-4-120-0001.png")
+    //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single.png")
+    //use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/rbc_single_oblong_2.png")
+    use img = new Image<Bgr, byte>("../../../../src/Tests_hough/images/two_rbc_1.png")
+   
+    let result = ImageAnalysis.doAnalysis img config
+
+    display mainWindow result.img
+    mainWindow.Root.Show()
+
+    app.Run() |> ignore
diff --git a/Parasitemia/Parasitemia/Types.fs b/Parasitemia/Parasitemia/Types.fs
new file mode 100644 (file)
index 0000000..308fc6b
--- /dev/null
@@ -0,0 +1,8 @@
+module Types
+
+open System
+
+type Ellipse = { cx: float; cy: float; a: float; b: float; alpha: float }
+    
+
+//type PointImg = { x: int; y: int }
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/Utils.fs b/Parasitemia/Parasitemia/Utils.fs
new file mode 100644 (file)
index 0000000..de16664
--- /dev/null
@@ -0,0 +1,19 @@
+module Utils
+
+open System.Diagnostics
+
+let roundInt = int << round
+
+let inline dprintfn fmt =
+    Printf.ksprintf System.Diagnostics.Debug.WriteLine fmt
+
+let mutable log : (string -> unit) =
+    fun m -> ()
+
+let logTime (m: string) (f: unit -> 'a) : 'a =
+    let sw = Stopwatch()
+    sw.Start()
+    let res = f ()
+    sw.Stop()
+    log <| sprintf "%A (time: %A ms)" m sw.ElapsedMilliseconds
+    res 
\ No newline at end of file
diff --git a/Parasitemia/Parasitemia/packages.config b/Parasitemia/Parasitemia/packages.config
new file mode 100644 (file)
index 0000000..1e68f7b
--- /dev/null
@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="utf-8"?>
+<packages>
+  <package id="Castle.Core" version="1.1.0" targetFramework="net46" />
+  <package id="Expression.Blend.Sdk" version="1.0.2" targetFramework="net46" />
+  <package id="FSharp.ViewModule.Core" version="0.9.9.1" targetFramework="net46" />
+  <package id="FsXaml.Wpf" version="0.9.9" targetFramework="net46" />
+  <package id="log4net" version="1.2.10" targetFramework="net46" />
+</packages>
\ No newline at end of file
diff --git a/Parasitemia/WPF/BitmapSourceConverter.cs b/Parasitemia/WPF/BitmapSourceConverter.cs
new file mode 100644 (file)
index 0000000..413a53b
--- /dev/null
@@ -0,0 +1,47 @@
+//----------------------------------------------------------------------------
+//  Copyright (C) 2004-2015 by EMGU Corporation. All rights reserved.       
+//----------------------------------------------------------------------------
+
+using System;
+using System.Runtime.InteropServices;
+using System.Windows;
+using System.Windows.Media;
+using System.Windows.Media.Imaging;
+using System.Windows.Shapes;
+using Emgu.CV;
+
+namespace Emgu.CV.WPF
+{
+   public static class BitmapSourceConvert
+   {
+      /// <summary>
+      /// Delete a GDI object
+      /// </summary>
+      /// <param name="o">The poniter to the GDI object to be deleted</param>
+      /// <returns></returns>
+      [DllImport("gdi32")]
+      private static extern int DeleteObject(IntPtr o);
+
+      /// <summary>
+      /// Convert an IImage to a WPF BitmapSource. The result can be used in the Set Property of Image.Source
+      /// </summary>
+      /// <param name="image">The Emgu CV Image</param>
+      /// <returns>The equivalent BitmapSource</returns>
+      public static BitmapSource ToBitmapSource(IImage image)
+      {
+         using (System.Drawing.Bitmap source = image.Bitmap)
+         {
+            IntPtr ptr = source.GetHbitmap(); //obtain the Hbitmap
+
+            BitmapSource bs = System.Windows.Interop.Imaging.CreateBitmapSourceFromHBitmap(
+                ptr,
+                IntPtr.Zero,
+                Int32Rect.Empty,
+                System.Windows.Media.Imaging.BitmapSizeOptions.FromEmptyOptions());
+
+            DeleteObject(ptr); //release the HBitmap
+            return bs;
+         }
+      }
+   }
+}
\ No newline at end of file
diff --git a/Parasitemia/WPF/Properties/AssemblyInfo.cs b/Parasitemia/WPF/Properties/AssemblyInfo.cs
new file mode 100644 (file)
index 0000000..d5a397f
--- /dev/null
@@ -0,0 +1,36 @@
+using System.Reflection;
+using System.Runtime.CompilerServices;
+using 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.
+[assembly: AssemblyTitle("WPF")]
+[assembly: AssemblyDescription("")]
+[assembly: AssemblyConfiguration("")]
+[assembly: AssemblyCompany("")]
+[assembly: AssemblyProduct("WPF")]
+[assembly: AssemblyCopyright("Copyright ©  2015")]
+[assembly: AssemblyTrademark("")]
+[assembly: AssemblyCulture("")]
+
+// 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.
+[assembly: ComVisible(false)]
+
+// The following GUID is for the ID of the typelib if this project is exposed to COM
+[assembly: Guid("314fd78e-870e-4794-bb16-ea4586f2abdb")]
+
+// 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:
+// [assembly: AssemblyVersion("1.0.*")]
+[assembly: AssemblyVersion("1.0.0.0")]
+[assembly: AssemblyFileVersion("1.0.0.0")]
diff --git a/Parasitemia/WPF/WPF.csproj b/Parasitemia/WPF/WPF.csproj
new file mode 100644 (file)
index 0000000..66a9cc9
--- /dev/null
@@ -0,0 +1,63 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="14.0" DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <Import Project="$(MSBuildExtensionsPath)\$(MSBuildToolsVersion)\Microsoft.Common.props" Condition="Exists('$(MSBuildExtensionsPath)\$(MSBuildToolsVersion)\Microsoft.Common.props')" />
+  <PropertyGroup>
+    <Configuration Condition=" '$(Configuration)' == '' ">Debug</Configuration>
+    <Platform Condition=" '$(Platform)' == '' ">AnyCPU</Platform>
+    <ProjectGuid>{314FD78E-870E-4794-BB16-EA4586F2ABDB}</ProjectGuid>
+    <OutputType>Library</OutputType>
+    <AppDesignerFolder>Properties</AppDesignerFolder>
+    <RootNamespace>WPF</RootNamespace>
+    <AssemblyName>WPF</AssemblyName>
+    <TargetFrameworkVersion>v4.0</TargetFrameworkVersion>
+    <FileAlignment>512</FileAlignment>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Debug|AnyCPU' ">
+    <DebugSymbols>true</DebugSymbols>
+    <DebugType>full</DebugType>
+    <Optimize>false</Optimize>
+    <OutputPath>bin\Debug\</OutputPath>
+    <DefineConstants>DEBUG;TRACE</DefineConstants>
+    <ErrorReport>prompt</ErrorReport>
+    <WarningLevel>4</WarningLevel>
+  </PropertyGroup>
+  <PropertyGroup Condition=" '$(Configuration)|$(Platform)' == 'Release|AnyCPU' ">
+    <DebugType>pdbonly</DebugType>
+    <Optimize>true</Optimize>
+    <OutputPath>bin\Release\</OutputPath>
+    <DefineConstants>TRACE</DefineConstants>
+    <ErrorReport>prompt</ErrorReport>
+    <WarningLevel>4</WarningLevel>
+  </PropertyGroup>
+  <ItemGroup>
+    <Reference Include="Emgu.CV">
+      <HintPath>..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.CV.dll</HintPath>
+    </Reference>
+    <Reference Include="Emgu.Util">
+      <HintPath>..\..\..\Emgu\emgucv-windows-universal 3.0.0.2157\bin\Emgu.Util.dll</HintPath>
+    </Reference>
+    <Reference Include="PresentationCore" />
+    <Reference Include="PresentationFramework" />
+    <Reference Include="System" />
+    <Reference Include="System.Core" />
+    <Reference Include="System.Drawing" />
+    <Reference Include="System.Xml.Linq" />
+    <Reference Include="System.Data.DataSetExtensions" />
+    <Reference Include="Microsoft.CSharp" />
+    <Reference Include="System.Data" />
+    <Reference Include="System.Xml" />
+    <Reference Include="WindowsBase" />
+  </ItemGroup>
+  <ItemGroup>
+    <Compile Include="BitmapSourceConverter.cs" />
+    <Compile Include="Properties\AssemblyInfo.cs" />
+  </ItemGroup>
+  <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
+  <!-- To modify your build process, add your task inside one of the targets below and uncomment it. 
+       Other similar extension points exist, see Microsoft.Common.targets.
+  <Target Name="BeforeBuild">
+  </Target>
+  <Target Name="AfterBuild">
+  </Target>
+  -->
+</Project>
\ No newline at end of file