Split the module 'ImgTools' in many modules.
authorGreg Burri <greg.burri@gmail.com>
Tue, 26 Jan 2016 13:55:21 +0000 (14:55 +0100)
committerGreg Burri <greg.burri@gmail.com>
Tue, 26 Jan 2016 13:55:21 +0000 (14:55 +0100)
17 files changed:
Parasitemia/ParasitemiaCore/Classifier.fs
Parasitemia/ParasitemiaCore/EEOver.fs
Parasitemia/ParasitemiaCore/Ellipse.fs
Parasitemia/ParasitemiaCore/Granulometry.fs
Parasitemia/ParasitemiaCore/ImgTools.fs [deleted file]
Parasitemia/ParasitemiaCore/ImgTools/Drawing.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/ImgTools/Edges.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/ImgTools/Histogram.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/ImgTools/IO.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/ImgTools/ImgTools.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/ImgTools/Otsu.fs [new file with mode: 0644]
Parasitemia/ParasitemiaCore/MainAnalysis.fs
Parasitemia/ParasitemiaCore/ParasitemiaCore.fsproj
Parasitemia/ParasitemiaCore/ParasitesMarker.fs
Parasitemia/ParasitemiaCore/UnitsOfMeasure.fs
Parasitemia/ParasitemiaUI/XAML/AboutWindow.xaml

index 615840c..f4b0844 100644 (file)
@@ -196,7 +196,6 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
                         then
                             parasiteArea <- parasiteArea + 1
 
-
                 let cellClass =
                     if float darkStainPixels > config.Parameters.maxDarkStainRatio * (float nbElement)
                     then
@@ -204,7 +203,7 @@ let findCells (ellipses: Ellipse list) (parasites: ParasitesMarker.Result) (img:
 
                     elif nucleusPixels.Count > 0 && parasiteArea >= minimumParasiteArea
                     then
-                        let infectionToRemove = ImgTools.connectedComponents parasites.parasite nucleusPixels
+                        let infectionToRemove = Morpho.connectedComponents parasites.parasite nucleusPixels
                         for p in infectionToRemove do
                             nucleusData.[p.Y, p.X, 0] <- 0uy
                         InfectedRBC
index eb8ac0b..6862dc5 100644 (file)
@@ -1,4 +1,5 @@
-module ParasitemiaCore.EEOver
+// Translation from https://github.com/chraibi/EEOver.
+module ParasitemiaCore.EEOver
 
 open System
 
@@ -508,7 +509,9 @@ let private biquadroots (p: float[]) (r: float[,]) =
 
     quad ()
 
-// Return a tuple (area, x intersections, y intersections)
+/// <summary>
+/// Return a tuple (area, x intersections, y intersections).
+/// </summary>
 let EEOverlapArea (e1: Types.Ellipse) (e2: Types.Ellipse) : (float32 * float32[] * float32[]) option =
     let h1, k1, a1, b1, phi_1 = float e1.Cx, float e1.Cy, float e1.A, float e1.B, float e1.Alpha
     let h2, k2, a2, b2, phi_2 = float e2.Cx, float e2.Cy, float e2.A, float e2.B, float e2.Alpha
index 1059e3b..4e48230 100644 (file)
@@ -14,142 +14,12 @@ open Config
 open MatchingEllipses
 open Const
 
-type private SearchExtremum = Minimum | Maximum
-
-let private goldenSectionSearch (f: float -> float) (nbIter: int) (xmin: float) (xmax: float) (searchExtremum: SearchExtremum) : (float * float) =
-    let gr = 1. / 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.
-    x, f x
-
-// Ellipse.A is always equal or greater than Ellipse.B.
-// Ellipse.Alpha is between 0 and Pi.
+/// <summary>
+/// Try to build an ellipse from three points and two tangents passing by the first and the second point.
+/// 'Ellipse.A' is always equal or greater than Ellipse.B.
+/// 'Ellipse.Alpha' is between 0 and Pi.
+/// </summary>
 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 = 10 // 3
-    let accuracy_extremum_search_2 = 10 // 4
-
-    // 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. + p1y ** 2.)
-    let theta1 = atan2 p1y p1x
-
-    let r2 = sqrt (p2x ** 2. + p2y ** 2.)
-    let theta2 = atan2 p2y p2x
-
-    let valid =
-        4. * 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. * sin (theta1 - theta2) ** 2. < 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. - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.)
-
-        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. 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. / 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.
-        let Wy = p1y + (p2y - p1y) / 2.
-
-        let Zx = p1x + (PTanx - p1x) / 2.
-        let Zy = p1y + (PTany - p1y) / 2.
-
-        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. + cy ** 2.)
-        let psi = atan2 cy cx
-
-        let rellipse theta =
-            sqrt (
-                rc ** 2. + (r1 ** 2. * r2 ** 2. * (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. * sin (theta1 - theta2) ** 2.) /
-                (sin (alpha1 - theta1) * sin (alpha2 - theta2) * (r1 * sin (theta - theta1) - r2 * sin (theta - theta2)) ** 2. - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.) ** 2. -
-                (2. * 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. - r1 * r2 * sin (alpha1 - theta) * sin (alpha2 - theta) * sin (theta1 - theta2) ** 2.))
-
-        // 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. (Math.PI / 2.) Maximum // Pi/2 and not pi because the period is Pi.
-        let r2eTheta, r2e = goldenSectionSearch rellipse accuracy_extremum_search_2 0. (Math.PI / 2.) 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.
-        then
-           alpha <- alpha + Math.PI
-
-        // Ride off the p3 referential.
-        let cx = cx + p3x
-        let cy = cy + p3y
-
-        Some (Types.Ellipse(float32 cx, float32 cy, float32 r1e, float32 r2e, float32 alpha))
-    else
-        None
-
-let ellipse2 (p1x: float) (p1y: float) (m1: float) (p2x: float) (p2y: float) (m2: float) (p3x: float) (p3y: float) : Types.Ellipse option =
     let p0 = pointFromTwoLines (Types.Line(float32 m1, float32 (p1y - m1 * p1x))) (Types.Line(float32 m2, float32(p2y - m2 * p2x)))
     let p0x, p0y = float p0.X, float p0.Y
 
@@ -259,7 +129,6 @@ let private areVectorsValid (p1x: float32) (p1y: float32) (p2x: float32) (p2y: f
         else
         Some (m1, m2)
 
-
 let find (edges: Matrix<byte>)
          (xGradient: Matrix<float32>)
          (yGradient: Matrix<float32>)
@@ -336,7 +205,7 @@ let find (edges: Matrix<byte>)
                         then
                             match areVectorsValid (float32 p1xf) (float32 p1yf) (float32 p2xf) (float32 p2yf) -xDirData.[p1.Y, p1.X] -yDirData.[p1.Y, p1.X] -xDirData.[p2.Y, p2.X] -yDirData.[p2.Y, p2.X] with
                             | Some (m1, m2) ->
-                                match ellipse2 p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf with
+                                match ellipse p1xf p1yf (float m1) p2xf p2yf (float m2) p3xf p3yf with
                                 | Some e when e.Cx > 0.f && e.Cx < w_f - 1.f && e.Cy > 0.f && e.Cy < h_f - 1.f &&
                                               e.A >= r1 - radiusTolerance && e.A <= r2 + radiusTolerance && e.B >= r1 - radiusTolerance && e.B <= r2 + radiusTolerance ->
                                      ellipses.Add e
index 0afe5cc..b325534 100644 (file)
@@ -71,7 +71,7 @@ let findRadiusByAreaClosing (img: Image<Gray, float32>) (radiusRange: int * int)
     let mutable maxDiff = 0.f
     let mutable max_r = r1
 
-    ImgTools.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (fun r diff ->
+    Morpho.areaCloseFWithFun imgCopy [ for r in r1 .. r2 -> Math.PI * float r ** 2. |> roundInt, r ] (fun r diff ->
         if r <> r1 && diff > maxDiff
         then
             maxDiff <- diff
diff --git a/Parasitemia/ParasitemiaCore/ImgTools.fs b/Parasitemia/ParasitemiaCore/ImgTools.fs
deleted file mode 100644 (file)
index cbb65b4..0000000
+++ /dev/null
@@ -1,1018 +0,0 @@
-module ParasitemiaCore.ImgTools
-
-open System
-open System.Drawing
-open System.Collections.Generic
-open System.Linq
-
-open Emgu.CV
-open Emgu.CV.Structure
-
-open Heap
-open Const
-open Types
-open Utils
-
-let normalize (img: Image<Gray, float32>) (upperLimit: float) : Image<Gray, float32> =
-    let min = ref [| 0.0 |]
-    let minLocation = ref <| [| Point() |]
-    let max = ref [| 0.0 |]
-    let maxLocation = ref <| [| Point() |]
-    img.MinMax(min, max, minLocation, maxLocation)
-    let normalized = (img - (!min).[0]) / ((!max).[0] - (!min).[0])
-    if upperLimit = 1.0
-    then normalized
-    else upperLimit * normalized
-
-let mergeChannels (img: Image<Bgr, float32>) (rgbWeights: float * float * float) : Image<Gray, float32> =
-    match rgbWeights with
-    | 1., 0., 0. -> img.[2]
-    | 0., 1., 0. -> img.[1]
-    | 0., 0., 1. -> img.[0]
-    | redFactor, greenFactor, blueFactor ->
-        let result = new Image<Gray, float32>(img.Size)
-        CvInvoke.AddWeighted(result, 1., img.[2], redFactor, 0., result)
-        CvInvoke.AddWeighted(result, 1., img.[1], greenFactor, 0., result)
-        CvInvoke.AddWeighted(result, 1., img.[0], blueFactor, 0., result)
-        result
-
-let mergeChannelsWithProjection (img: Image<Bgr, float32>) (v1r: float32, v1g: float32, v1b: float32) (v2r: float32, v2g: float32, v2b: float32) (upperLimit: float)  : Image<Gray, float32> =
-    let vr, vg, vb = v2r - v1r, v2g - v1g, v2b - v1b
-    let vMagnitude = sqrt (vr ** 2.f + vg ** 2.f + vb ** 2.f)
-    let project (r: float32) (g: float32) (b: float32) = ((r - v1r) * vr + (g - v1g) * vg + (b - v1b) * vb) / vMagnitude
-    let result = new Image<Gray, float32>(img.Size)
-    // TODO: Essayer en bindant Data pour gagner du temps
-    for i in 0 .. img.Height - 1 do
-        for j in 0 .. img.Width - 1 do
-            result.Data.[i, j, 0] <- project img.Data.[i, j, 2] img.Data.[i, j, 1] img.Data.[i, j, 0]
-    normalize result upperLimit
-
-// Normalize image values between 0uy and 255uy.
-let normalizeAndConvert (img: Image<Gray, 'TDepth>) : Image<Gray, byte> =
-    (normalize (img.Convert<Gray, float32>()) 255.).Convert<Gray, byte>()
-
-let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
-    img.Save(filepath)
-
-let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
-    use img = new Image<Gray, 'TDeph>(mat.Size)
-    mat.CopyTo(img)
-    saveImg img filepath
-
-type Histogram = { data: int[]; total: int; sum: int; min: float32; max: float32 }
-
-let histogramImg (img: Image<Gray, float32>) (nbSamples: int) : Histogram =
-    let imgData = img.Data
-
-    let min, max =
-        let min = ref [| 0.0 |]
-        let minLocation = ref <| [| Point() |]
-        let max = ref [| 0.0 |]
-        let maxLocation = ref <| [| Point() |]
-        img.MinMax(min, max, minLocation, maxLocation)
-        float32 (!min).[0], float32 (!max).[0]
-
-    let inline bin (x: float32) : int =
-        let p = int ((x - min) / (max - min) * float32 nbSamples)
-        if p >= nbSamples then nbSamples - 1 else p
-
-    let data = Array.zeroCreate nbSamples
-
-    for i in 0 .. img.Height - 1 do
-        for j in 0 .. img.Width - 1 do
-            let p = bin imgData.[i, j, 0]
-            data.[p] <- data.[p] + 1
-
-    { data = data; total = img.Height * img.Width; sum = Array.sum data; min = min; max = max }
-
-let histogramMat (mat: Matrix<float32>) (nbSamples: int) : Histogram =
-    let matData = mat.Data
-
-    let min, max =
-        let min = ref 0.0
-        let minLocation = ref <| Point()
-        let max = ref 0.0
-        let maxLocation = ref <| Point()
-        mat.MinMax(min, max, minLocation, maxLocation)
-        float32 !min, float32 !max
-
-    let inline bin (x: float32) : int =
-        let p = int ((x - min) / (max - min) * float32 nbSamples)
-        if p >= nbSamples then nbSamples - 1 else p
-
-    let data = Array.zeroCreate nbSamples
-
-    for i in 0 .. mat.Height - 1 do
-        for j in 0 .. mat.Width - 1 do
-            let p = bin matData.[i, j]
-            data.[p] <- data.[p] + 1
-
-    { data = data; total = mat.Height * mat.Width; sum = Array.sum data; min = min; max = max }
-
-let histogram (values: float32 seq) (nbSamples: int) : Histogram =
-    let mutable min = Single.MaxValue
-    let mutable max = Single.MinValue
-    let mutable n = 0
-
-    for v in values do
-        n <- n + 1
-        if v < min then min <- v
-        if v > max then max <- v
-
-    let inline bin (x: float32) : int =
-        let p = int ((x - min) / (max - min) * float32 nbSamples)
-        if p >= nbSamples then nbSamples - 1 else p
-
-    let data = Array.zeroCreate nbSamples
-
-    for v in values do
-        let p = bin v
-        data.[p] <- data.[p] + 1
-
-    { data = data; total = n; sum = Array.sum data; min = min; max = max }
-
-let otsu (hist: Histogram) : float32 * float32 * float32 =
-    let mutable sumB = 0
-    let mutable wB = 0
-    let mutable maximum = 0.0
-    let mutable level = 0
-    let sum = hist.data |> Array.mapi (fun i v -> i * v |> float) |> Array.sum
-
-    for i in 0 .. hist.data.Length - 1 do
-        wB <- wB + hist.data.[i]
-        if wB <> 0
-        then
-            let wF = hist.total - wB
-            if wF <> 0
-            then
-                sumB <- sumB + i * hist.data.[i]
-                let mB = (float sumB) / (float wB)
-                let mF = (sum - float sumB) / (float wF)
-                let between = (float wB) * (float wF) * (mB - mF) ** 2.;
-                if between >= maximum
-                then
-                    level <- i
-                    maximum <- between
-
-    let mean1 =
-        let mutable sum = 0
-        let mutable nb = 0
-        for i in 0 .. level - 1 do
-            sum <- sum + i * hist.data.[i]
-            nb <- nb + hist.data.[i]
-        (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2)
-
-    let mean2 =
-        let mutable sum = 0
-        let mutable nb = 0
-        for i in level + 1 .. hist.data.Length - 1 do
-            sum <- sum + i * hist.data.[i]
-            nb <- nb + hist.data.[i]
-        (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2)
-
-    let toFloat l =
-        float32 l / float32 hist.data.Length * (hist.max - hist.min) + hist.min
-
-    toFloat level, toFloat mean1, toFloat mean2
-
-/// <summary>
-/// Remove M-adjacent pixels. It may be used after thinning.
-/// </summary>
-let suppressMAdjacency (img: Matrix<byte>) =
-    let w = img.Width
-    let h = img.Height
-    for i in 1 .. h - 2 do
-        for j in 1 .. w - 2 do
-            if img.[i, j] > 0uy && img.Data.[i + 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i - 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i - 1, j - 1] = 0uy)
-            then
-                img.[i, j] <- 0uy
-    for i in 1 .. h - 2 do
-        for j in 1 .. w - 2 do
-            if img.[i, j] > 0uy && img.Data.[i - 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i + 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i + 1, j - 1] = 0uy)
-            then
-                img.[i, j] <- 0uy
-
-/// <summary>
-/// Find edges of an image by using the Canny approach.
-/// The thresholds are automatically defined with otsu on gradient magnitudes.
-/// </summary>
-/// <param name="img"></param>
-let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Matrix<float32> * Matrix<float32> =
-    let w = img.Width
-    let h = img.Height
-
-    use sobelKernel =
-        new Matrix<float32>(array2D [[ 1.0f; 0.0f; -1.0f ]
-                                     [ 2.0f; 0.0f; -2.0f ]
-                                     [ 1.0f; 0.0f; -1.0f ]])
-
-    let xGradient = new Matrix<float32>(img.Size)
-    let yGradient = new Matrix<float32>(img.Size)
-    CvInvoke.Filter2D(img, xGradient, sobelKernel, Point(1, 1))
-    CvInvoke.Filter2D(img, yGradient, sobelKernel.Transpose(), Point(1, 1))
-
-    use magnitudes = new Matrix<float32>(xGradient.Size)
-    use angles = new Matrix<float32>(xGradient.Size)
-    CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles.
-
-    let thresholdHigh, thresholdLow =
-        let sensibilityHigh = 0.1f
-        let sensibilityLow = 0.0f
-        let threshold, _, _ = otsu (histogramMat magnitudes 300)
-        threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
-
-    // Non-maximum suppression.
-    use nms = new Matrix<byte>(xGradient.Size)
-
-    let nmsData = nms.Data
-    let anglesData = angles.Data
-    let magnitudesData = magnitudes.Data
-    let xGradientData = xGradient.Data
-    let yGradientData = yGradient.Data
-
-    for i in 0 .. h - 1 do
-        nmsData.[i, 0] <- 0uy
-        nmsData.[i, w - 1] <- 0uy
-
-    for j in 0 .. w - 1 do
-        nmsData.[0, j] <- 0uy
-        nmsData.[h - 1, j] <- 0uy
-
-    for i in 1 .. h - 2 do
-        for j in 1 .. w - 2 do
-            let vx = xGradientData.[i, j]
-            let vy = yGradientData.[i, j]
-            if vx <> 0.f || vy <> 0.f
-            then
-                let angle = anglesData.[i, j]
-
-                let vx', vy' = abs vx, abs vy
-                let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
-                let ratio1 = 1.f - ratio2
-
-                let mNeigbors (sign: int) : float32 =
-                    if angle < PI / 4.f
-                    then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
-                    elif angle < PI / 2.f
-                    then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j]
-                    elif angle < 3.f * PI / 4.f
-                    then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign]
-                    elif angle < PI
-                    then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign]
-                    elif angle < 5.f * PI / 4.f
-                    then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign]
-                    elif angle < 3.f * PI / 2.f
-                    then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j]
-                    elif angle < 7.f * PI / 4.f
-                    then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign]
-                    else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]
-
-                let m = magnitudesData.[i, j]
-                if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
-                then
-                    nmsData.[i, j] <- 1uy
-
-    // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
-
-    let edges = new Matrix<byte>(xGradient.Size)
-    let edgesData = edges.Data
-
-    // Hysteresis thresholding.
-    let toVisit = Stack<Point>()
-    for i in 0 .. h - 1 do
-        for j in 0 .. w - 1 do
-            if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh
-            then
-                nmsData.[i, j] <- 0uy
-                toVisit.Push(Point(j, i))
-                while toVisit.Count > 0 do
-                    let p = toVisit.Pop()
-                    edgesData.[p.Y, p.X] <- 1uy
-                    for i' in -1 .. 1  do
-                        for j' in -1 .. 1 do
-                            if i' <> 0 || j' <> 0
-                            then
-                                let ni = p.Y + i'
-                                let nj = p.X + j'
-                                if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
-                                then
-                                    nmsData.[ni, nj] <- 0uy
-                                    toVisit.Push(Point(nj, ni))
-
-    edges, xGradient, yGradient
-
-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)
-
-let drawPoints (img: Image<Gray, 'TDepth>) (points: Points) (intensity: 'TDepth) =
-    for p in points do
-        img.Data.[p.Y, p.X, 0] <- intensity
-
-type ExtremumType =
-    | Maxima = 1
-    | Minima = 2
-
-let findExtremum (img: Image<Gray, 'TDepth>) (extremumType: ExtremumType) : IEnumerable<Points> =
-    let w = img.Width
-    let h = img.Height
-    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
-
-    let imgData = img.Data
-    let suppress: bool[,] = Array2D.zeroCreate h w
-
-    let result = List<List<Point>>()
-
-    let flood (start: Point) : List<List<Point>> =
-        let sameLevelToCheck = Stack<Point>()
-        let betterLevelToCheck = Stack<Point>()
-        betterLevelToCheck.Push(start)
-
-        let result' = List<List<Point>>()
-
-        while betterLevelToCheck.Count > 0 do
-            let p = betterLevelToCheck.Pop()
-            if not suppress.[p.Y, p.X]
-            then
-                suppress.[p.Y, p.X] <- true
-                sameLevelToCheck.Push(p)
-                let current = List<Point>()
-
-                let mutable betterExists = false
-
-                while sameLevelToCheck.Count > 0 do
-                    let p' = sameLevelToCheck.Pop()
-                    let currentLevel = imgData.[p'.Y, p'.X, 0]
-                    current.Add(p') |> ignore
-                    for i, j in se do
-                        let ni = i + p'.Y
-                        let nj = j + p'.X
-                        if ni >= 0 && ni < h && nj >= 0 && nj < w
-                        then
-                            let level = imgData.[ni, nj, 0]
-                            let notSuppressed = not suppress.[ni, nj]
-
-                            if level = currentLevel && notSuppressed
-                            then
-                                suppress.[ni, nj] <- true
-                                sameLevelToCheck.Push(Point(nj, ni))
-                            elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel
-                            then
-                                betterExists <- true
-                                if notSuppressed
-                                then
-                                    betterLevelToCheck.Push(Point(nj, ni))
-
-                if not betterExists
-                then
-                    result'.Add(current)
-        result'
-
-    for i in 0 .. h - 1 do
-        for j in 0 .. w - 1 do
-            let maxima = flood (Point(j, i))
-            if maxima.Count > 0
-            then
-                result.AddRange(maxima)
-
-    result.Select(fun l -> Points(l))
-
-let findMaxima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
-    findExtremum img ExtremumType.Maxima
-
-let findMinima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
-    findExtremum img ExtremumType.Minima
-
-type PriorityQueue () =
-    let size = 256
-    let q: Points[] = Array.init size (fun i -> Points())
-    let mutable highest = -1 // Value of the first elements of 'q'.
-    let mutable lowest = size
-
-    member this.NextMax () : byte * Point =
-        if this.IsEmpty
-        then
-            invalidOp "Queue is empty"
-        else
-            let l = q.[highest]
-            let next = l.First()
-            l.Remove(next) |> ignore
-            let value = byte highest
-
-            if l.Count = 0
-            then
-                highest <- highest - 1
-                while highest > lowest && q.[highest].Count = 0 do
-                    highest <- highest - 1
-                if highest = lowest
-                then
-                    highest <- -1
-                    lowest <- size
-
-            value, next
-
-    member this.NextMin () : byte * Point =
-        if this.IsEmpty
-        then
-            invalidOp "Queue is empty"
-        else
-            let l = q.[lowest + 1]
-            let next = l.First()
-            l.Remove(next) |> ignore
-            let value = byte (lowest + 1)
-
-            if l.Count = 0
-            then
-                lowest <- lowest + 1
-                while lowest < highest && q.[lowest + 1].Count = 0 do
-                    lowest <- lowest + 1
-                if highest = lowest
-                then
-                    highest <- -1
-                    lowest <- size
-
-            value, next
-
-    member this.Max =
-        highest |> byte
-
-    member this.Min =
-        lowest + 1 |> byte
-
-    member this.Add (value: byte) (p: Point) =
-        let vi = int value
-
-        if vi > highest
-        then
-            highest <- vi
-        if vi <= lowest
-        then
-            lowest <- vi - 1
-
-        q.[vi].Add(p) |> ignore
-
-    member this.Remove (value: byte) (p: Point) =
-        let vi = int value
-        if q.[vi].Remove(p) && q.[vi].Count = 0
-        then
-            if vi = highest
-            then
-                highest <- highest - 1
-                while highest > lowest && q.[highest].Count = 0 do
-                    highest <- highest - 1
-            elif vi - 1 = lowest
-            then
-                lowest <- lowest + 1
-                while lowest < highest && q.[lowest + 1].Count = 0 do
-                    lowest <- lowest + 1
-
-            if highest = lowest // The queue is now empty.
-            then
-                highest <- -1
-                lowest <- size
-
-    member this.IsEmpty =
-        highest = -1
-
-    member this.Clear () =
-        while highest > lowest  do
-            q.[highest].Clear()
-            highest <- highest - 1
-        highest <- -1
-        lowest <- size
-
-type private AreaState =
-    | Removed = 1
-    | Unprocessed = 2
-    | Validated = 3
-
-type private AreaOperation =
-    | Opening = 1
-    | Closing = 2
-
-[<AllowNullLiteral>]
-type private Area (elements: Points) =
-    member this.Elements = elements
-    member val Intensity = None with get, set
-    member val State = AreaState.Unprocessed with get, set
-
-let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperation) =
-    let w = img.Width
-    let h = img.Height
-    let imgData = img.Data
-    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
-
-    let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
-
-    let pixels: Area[,] = Array2D.create h w null
-    for m in areas do
-        for e in m.Elements do
-            pixels.[e.Y, e.X] <- m
-
-    let queue = PriorityQueue()
-
-    let addEdgeToQueue (elements: Points) =
-        for p in elements do
-            for i, j in se do
-                let ni = i + p.Y
-                let nj = j + p.X
-                let p' = Point(nj, ni)
-                if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p'))
-                then
-                    queue.Add (imgData.[ni, nj, 0]) p'
-
-    // Reverse order is quicker.
-    for i in areas.Count - 1 .. -1 .. 0 do
-        let m = areas.[i]
-        if m.Elements.Count <= area && m.State <> AreaState.Removed
-        then
-            queue.Clear()
-            addEdgeToQueue m.Elements
-
-            let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min
-            let nextElements = Points()
-
-            let mutable stop = false
-            while not stop do
-                let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
-                let mutable merged = false
-
-                if intensity' = intensity // The intensity doesn't change.
-                then
-                    if m.Elements.Count + nextElements.Count + 1 > area
-                    then
-                        m.State <- AreaState.Validated
-                        m.Intensity <- Some intensity
-                        stop <- true
-                    else
-                        nextElements.Add(p) |> ignore
-
-                elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity
-                then
-                    m.Elements.UnionWith(nextElements)
-                    for e in nextElements do
-                        pixels.[e.Y, e.X] <- m
-
-                    if m.Elements.Count = area
-                    then
-                        m.State <- AreaState.Validated
-                        m.Intensity <- Some (intensity')
-                        stop <- true
-                    else
-                        intensity <- intensity'
-                        nextElements.Clear()
-                        nextElements.Add(p) |> ignore
-
-                else
-                    match pixels.[p.Y, p.X] with
-                    | null -> ()
-                    | m' ->
-                        if m'.Elements.Count + m.Elements.Count <= area
-                        then
-                            m'.State <- AreaState.Removed
-                            for e in m'.Elements do
-                                pixels.[e.Y, e.X] <- m
-                                queue.Remove imgData.[e.Y, e.X, 0] e
-                            addEdgeToQueue m'.Elements
-                            m.Elements.UnionWith(m'.Elements)
-                            let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min
-                            if intensityMax <> intensity
-                            then
-                                intensity <- intensityMax
-                                nextElements.Clear()
-                            merged <- true
-
-                    if not merged
-                    then
-                        m.State <- AreaState.Validated
-                        m.Intensity <- Some (intensity)
-                        stop <- true
-
-                if not stop && not merged
-                then
-                    for i, j in se do
-                        let ni = i + p.Y
-                        let nj = j + p.X
-                        let p' = Point(nj, ni)
-                        if ni < 0 || ni >= h || nj < 0 || nj >= w
-                        then
-                            m.State <- AreaState.Validated
-                            m.Intensity <- Some (intensity)
-                            stop <- true
-                        elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
-                        then
-                            queue.Add (imgData.[ni, nj, 0]) p'
-
-                if queue.IsEmpty
-                then
-                    if m.Elements.Count + nextElements.Count <= area
-                    then
-                        m.State <- AreaState.Validated
-                        m.Intensity <- Some intensity'
-                        m.Elements.UnionWith(nextElements)
-                    stop <- true
-
-    for m in areas do
-        if m.State = AreaState.Validated
-        then
-            match m.Intensity with
-            | Some i ->
-                for p in m.Elements do
-                    imgData.[p.Y, p.X, 0] <- i
-            | _ -> ()
-    ()
-
-/// <summary>
-/// Area opening on byte image.
-/// </summary>
-let areaOpen (img: Image<Gray, byte>) (area: int) =
-    areaOperation img area AreaOperation.Opening
-
-/// <summary>
-/// Area closing on byte image.
-/// </summary>
-let areaClose (img: Image<Gray, byte>) (area: int) =
-    areaOperation img area AreaOperation.Closing
-
-// A simpler algorithm than 'areaOpen' on byte image but slower.
-let areaOpen2 (img: Image<Gray, byte>) (area: int) =
-    let w = img.Width
-    let h = img.Height
-    let imgData = img.Data
-    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
-
-    let histogram = Array.zeroCreate 256
-    for i in 0 .. h - 1 do
-        for j in 0 .. w - 1 do
-            let v = imgData.[i, j, 0] |> int
-            histogram.[v] <- histogram.[v] + 1
-
-    let flooded : bool[,] = Array2D.zeroCreate h w
-
-    let pointsChecked = HashSet<Point>()
-    let pointsToCheck = Stack<Point>()
-
-    for level in 255 .. -1 .. 0 do
-        let mutable n = histogram.[level]
-        if n > 0
-        then
-            for i in 0 .. h - 1 do
-                for j in 0 .. w - 1 do
-                    if not flooded.[i, j] && imgData.[i, j, 0] = byte level
-                    then
-                        let mutable maxNeighborValue = 0uy
-                        pointsChecked.Clear()
-                        pointsToCheck.Clear()
-                        pointsToCheck.Push(Point(j, i))
-
-                        while pointsToCheck.Count > 0 do
-                            let next = pointsToCheck.Pop()
-                            pointsChecked.Add(next) |> ignore
-                            flooded.[next.Y, next.X] <- true
-
-                            for nx, ny in se do
-                                let p = Point(next.X + nx, next.Y + ny)
-                                if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
-                                then
-                                    let v = imgData.[p.Y, p.X, 0]
-                                    if v = byte level
-                                    then
-                                        if not (pointsChecked.Contains(p))
-                                        then
-                                            pointsToCheck.Push(p)
-                                    elif v > maxNeighborValue
-                                    then
-                                        maxNeighborValue <- v
-
-                        if int maxNeighborValue < level && pointsChecked.Count <= area
-                        then
-                            for p in pointsChecked do
-                                imgData.[p.Y, p.X, 0] <- maxNeighborValue
-
-[<AllowNullLiteral>]
-type Island (cmp: IComparer<float32>) =
-    member val Shore = Heap.Heap<float32, Point>(cmp) with get
-    member val Level = 0.f with get, set
-    member val Surface = 0 with get, set
-    member this.IsInfinite = this.Surface = Int32.MaxValue
-
-let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) =
-    let w = img.Width
-    let h = img.Height
-    let earth = img.Data
-    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
-
-    let comparer = if op = AreaOperation.Opening
-                   then { new IComparer<float32> with member this.Compare(v1, v2) = v1.CompareTo(v2) }
-                   else { new IComparer<float32> with member this.Compare(v1, v2) = v2.CompareTo(v1) }
-
-    let ownership: Island[,] = Array2D.create h w null
-
-    // Initialize islands with their shore.
-    let islands = List<Island>()
-    let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima
-    for e in extremum do
-        let island =
-            let p = e.First()
-            Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count)
-        islands.Add(island)
-        let shorePoints = Points()
-        for p in e do
-            ownership.[p.Y, p.X] <- island
-            for i, j in se do
-                let ni = i + p.Y
-                let nj = j + p.X
-                let neighbor = Point(nj, ni)
-                if ni >= 0 && ni < h && nj >= 0 && nj < w && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor))
-                then
-                    shorePoints.Add(neighbor) |> ignore
-                    island.Shore.Add earth.[ni, nj, 0] neighbor
-
-    for area, obj in areas do
-        for island in islands do
-            let mutable stop = island.Shore.IsEmpty
-
-            // 'true' if 'p' is owned or adjacent to 'island'.
-            let inline ownedOrAdjacent (p: Point) : bool =
-                ownership.[p.Y, p.X] = island ||
-                (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) ||
-                (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) ||
-                (p.X > 0 && ownership.[p.Y, p.X - 1] = island) ||
-                (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island)
-
-            while not stop && island.Surface < area do
-                let level, next = island.Shore.Max
-                let other = ownership.[next.Y, next.X]
-                if other = island // During merging, some points on the shore may be owned by the island itself -> ignored.
-                then
-                    island.Shore.RemoveNext ()
-                else
-                    if not <| Object.ReferenceEquals(other, null)
-                    then // We touching another island.
-                        if island.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area || comparer.Compare(island.Level, other.Level) < 0
-                        then
-                            stop <- true
-                        else // We can merge 'other' into 'surface'.
-                            island.Surface <- island.Surface + other.Surface
-                            island.Level <- other.Level
-                            // island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then other.Level else island.Level
-                            for l, p in other.Shore do
-                                let mutable currentY = p.Y + 1
-                                while currentY < h && ownership.[currentY, p.X] = other do
-                                    ownership.[currentY, p.X] <- island
-                                    currentY <- currentY + 1
-                                island.Shore.Add l p
-                            other.Shore.Clear()
-
-                    elif comparer.Compare(level, island.Level) > 0
-                    then
-                        stop <- true
-                    else
-                        island.Shore.RemoveNext ()
-                        for i, j in se do
-                            let ni = i + next.Y
-                            let nj = j + next.X
-                            if ni < 0 || ni >= h || nj < 0 || nj >= w
-                            then
-                                island.Surface <- Int32.MaxValue
-                                stop <- true
-                            else
-                                let neighbor = Point(nj, ni)
-                                if not <| ownedOrAdjacent neighbor
-                                then
-                                    island.Shore.Add earth.[ni, nj, 0] neighbor
-                        if not stop
-                        then
-                            ownership.[next.Y, next.X] <- island
-                            island.Level <- level
-                            island.Surface <- island.Surface + 1
-
-        let mutable diff = 0.f
-
-        for i in 0 .. h - 1 do
-            for j in 0 .. w - 1 do
-                match ownership.[i, j] with
-                | null -> ()
-                | island ->
-                    let l = island.Level
-                    diff <- diff + l - earth.[i, j, 0]
-                    earth.[i, j, 0] <- l
-
-        match f with
-        | Some f' -> f' obj diff
-        | _ -> ()
-    ()
-
-/// <summary>
-/// Area opening on float image.
-/// </summary>
-let areaOpenF (img: Image<Gray, float32>) (area: int) =
-    areaOperationF img [ area, () ] None AreaOperation.Opening
-
-/// <summary>
-/// Area closing on float image.
-/// </summary>
-let areaCloseF (img: Image<Gray, float32>) (area: int) =
-    areaOperationF img [ area, () ] None AreaOperation.Closing
-
-/// <summary>
-/// Area closing on float image with different areas. Given areas must be sorted increasingly.
-/// For each area the function 'f' is called with the associated area value of type 'a and the volume difference
-/// Between the previous and the current closing.
-/// </summary>
-let areaOpenFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
-    areaOperationF img areas (Some f) AreaOperation.Opening
-
-/// <summary>
-/// Same as 'areaOpenFWithFun' for closing operation.
-/// </summary>
-let areaCloseFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
-    areaOperationF img areas (Some f) AreaOperation.Closing
-
-/// <summary>
-/// Zhang and Suen thinning algorithm.
-/// Modify 'mat' in place.
-/// </summary>
-let thin (mat: Matrix<byte>) =
-    let w = mat.Width
-    let h = mat.Height
-    let mutable data1 = mat.Data
-    let mutable data2 = Array2D.copy data1
-
-    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 p2 = if i = 0 then 0uy else data1.[i-1, j]
-                    let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
-                    let p4 = if j = w-1 then 0uy else data1.[i, j+1]
-                    let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
-                    let p6 = if i = h-1 then 0uy else data1.[i+1, j]
-                    let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
-                    let p8 = if j = 0 then 0uy else data1.[i, j-1]
-                    let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
-
-                    let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
-                    if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
-                        (if p2 = 0uy && p3 = 1uy then 1 else 0) +
-                        (if p3 = 0uy && p4 = 1uy then 1 else 0) +
-                        (if p4 = 0uy && p5 = 1uy then 1 else 0) +
-                        (if p5 = 0uy && p6 = 1uy then 1 else 0) +
-                        (if p6 = 0uy && p7 = 1uy then 1 else 0) +
-                        (if p7 = 0uy && p8 = 1uy then 1 else 0) +
-                        (if p8 = 0uy && p9 = 1uy then 1 else 0) +
-                        (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
-                        if oddIteration
-                        then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
-                        else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
-                    then
-                        data2.[i, j] <- 0uy
-                        pixelChanged <- true
-                else
-                    data2.[i, j] <- 0uy
-
-        oddIteration <- not oddIteration
-        let tmp = data1
-        data1 <- data2
-        data2 <- tmp
-
-/// <summary>
-/// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
-/// Modify 'mat' in place.
-/// </summary>
-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
-
-    use 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<Point>()
-                let neighborsToCheck = Stack<Point>()
-                neighborsToCheck.Push(Point(j, i))
-                data'.[i, j] <- 0uy
-
-                while neighborsToCheck.Count > 0 do
-                    let n = neighborsToCheck.Pop()
-                    neighborhood.Add(n)
-                    for (ni, nj) in neighbors do
-                        let pi = n.Y + ni
-                        let pj = n.X + nj
-                        if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
-                        then
-                            neighborsToCheck.Push(Point(pj, pi))
-                            data'.[pi, pj] <- 0uy
-                if neighborhood.Count <= areaSize
-                then
-                    for n in neighborhood do
-                        data.[n.Y, n.X] <- 0uy
-
-let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Points =
-    let w = img.Width
-    let h = img.Height
-
-    let pointChecked = Points()
-    let pointToCheck = Stack<Point>(startPoints);
-
-    let data = img.Data
-
-    while pointToCheck.Count > 0 do
-        let next = pointToCheck.Pop()
-        pointChecked.Add(next) |> ignore
-        for ny in -1 .. 1 do
-            for nx in -1 .. 1 do
-                if ny <> 0 && nx <> 0
-                then
-                    let p = Point(next.X + nx, next.Y + ny)
-                    if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
-                    then
-                        pointToCheck.Push(p)
-
-    pointChecked
-
-let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
-    img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
-
-let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
-    img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
-
-let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Ellipse) (color: 'TColor) (alpha: float) =
-    if alpha >= 1.0
-    then
-        img.Draw(Emgu.CV.Structure.Ellipse(PointF(e.Cx, e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
-    else
-        let windowPosX = e.Cx - e.A - 5.f
-        let gapX = windowPosX - (float32 (int windowPosX))
-
-        let windowPosY = e.Cy - e.A - 5.f
-        let gapY = windowPosY - (float32 (int windowPosY))
-
-        let roi = Rectangle(int windowPosX, int windowPosY, 2.f * (e.A + 5.f) |> int, 2.f * (e.A + 5.f) |> int)
-
-        img.ROI <- roi
-        if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
-        then
-            use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
-            i.Draw(Emgu.CV.Structure.Ellipse(PointF(e.A + 5.f + gapX, e.A + 5.f + gapY), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
-            CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
-        img.ROI <- Rectangle.Empty
-
-let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Ellipse list) (color: 'TColor) (alpha: float) =
-    List.iter (fun e -> drawEllipse img e color alpha) ellipses
-
-let rngCell =  System.Random()
-let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Cell) =
-    if drawCellContent
-    then
-        let colorB = rngCell.Next(20, 70)
-        let colorG = rngCell.Next(20, 70)
-        let colorR = rngCell.Next(20, 70)
-
-        for y in 0 .. c.elements.Height - 1 do
-            for x in 0 .. c.elements.Width - 1 do
-                if c.elements.[y, x] > 0uy
-                then
-                    let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2
-                    let b = img.Data.[y + dy, x + dx, 0] |> int
-                    let g = img.Data.[y + dy, x + dx, 1] |> int
-                    let r = img.Data.[y + dy, x + dx, 2] |> int
-                    img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
-                    img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
-                    img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
-
-    let crossColor, crossColor2 =
-        match c.cellClass with
-        | HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
-        | InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
-        | Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
-
-    drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2
-    drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2
-
-    drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1
-    drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
-
-
-let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Cell list) =
-    List.iter (fun c -> drawCell img drawCellContent c) cells
\ No newline at end of file
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Drawing.fs b/Parasitemia/ParasitemiaCore/ImgTools/Drawing.fs
new file mode 100644 (file)
index 0000000..363acda
--- /dev/null
@@ -0,0 +1,79 @@
+module ParasitemiaCore.Drawing
+
+open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+open Const
+open Types
+
+let drawPoints (img: Image<Gray, 'TDepth>) (points: Points) (intensity: 'TDepth) =
+    for p in points do
+        img.Data.[p.Y, p.X, 0] <- intensity
+
+let drawLine (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: int) (y0: int) (x1: int) (y1: int) (thickness: int) =
+    img.Draw(LineSegment2D(Point(x0, y0), Point(x1, y1)), color, thickness);
+
+let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: float) (x1: float) (y1: float) (thickness: int) =
+    img.Draw(LineSegment2DF(PointF(float32 x0, float32 y0), PointF(float32 x1, float32 y1)), color, thickness, CvEnum.LineType.AntiAlias);
+
+let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Ellipse) (color: 'TColor) (alpha: float) =
+    if alpha >= 1.0
+    then
+        img.Draw(Emgu.CV.Structure.Ellipse(PointF(e.Cx, e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
+    else
+        let windowPosX = e.Cx - e.A - 5.f
+        let gapX = windowPosX - (float32 (int windowPosX))
+
+        let windowPosY = e.Cy - e.A - 5.f
+        let gapY = windowPosY - (float32 (int windowPosY))
+
+        let roi = Rectangle(int windowPosX, int windowPosY, 2.f * (e.A + 5.f) |> int, 2.f * (e.A + 5.f) |> int)
+
+        img.ROI <- roi
+        if roi = img.ROI // We do not display ellipses touching the edges (FIXME)
+        then
+            use i = new Image<'TColor, 'TDepth>(img.ROI.Size)
+            i.Draw(Emgu.CV.Structure.Ellipse(PointF(e.A + 5.f + gapX, e.A + 5.f + gapY), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
+            CvInvoke.AddWeighted(img, 1.0, i, alpha, 0.0, img)
+        img.ROI <- Rectangle.Empty
+
+let drawEllipses (img: Image<'TColor, 'TDepth>) (ellipses: Ellipse list) (color: 'TColor) (alpha: float) =
+    List.iter (fun e -> drawEllipse img e color alpha) ellipses
+
+let rngCell =  System.Random()
+let drawCell (img: Image<Bgr, byte>) (drawCellContent: bool) (c: Cell) =
+    if drawCellContent
+    then
+        let colorB = rngCell.Next(20, 70)
+        let colorG = rngCell.Next(20, 70)
+        let colorR = rngCell.Next(20, 70)
+
+        for y in 0 .. c.elements.Height - 1 do
+            for x in 0 .. c.elements.Width - 1 do
+                if c.elements.[y, x] > 0uy
+                then
+                    let dx, dy = c.center.X - c.elements.Width / 2, c.center.Y - c.elements.Height / 2
+                    let b = img.Data.[y + dy, x + dx, 0] |> int
+                    let g = img.Data.[y + dy, x + dx, 1] |> int
+                    let r = img.Data.[y + dy, x + dx, 2] |> int
+                    img.Data.[y + dy, x + dx, 0] <- if b + colorB > 255 then 255uy else byte (b + colorB)
+                    img.Data.[y + dy, x + dx, 1] <- if g + colorG > 255 then 255uy else byte (g + colorG)
+                    img.Data.[y + dy, x + dx, 2] <- if r + colorR > 255 then 255uy else byte (r + colorR)
+
+    let crossColor, crossColor2 =
+        match c.cellClass with
+        | HealthyRBC -> Bgr(255., 0., 0.), Bgr(255., 255., 255.)
+        | InfectedRBC -> Bgr(0., 0., 255.), Bgr(120., 120., 255.)
+        | Peculiar -> Bgr(0., 0., 0.), Bgr(80., 80., 80.)
+
+    drawLine img crossColor2 (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 2
+    drawLine img crossColor2 c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 2
+
+    drawLine img crossColor (c.center.X - 3) c.center.Y (c.center.X + 3) c.center.Y 1
+    drawLine img crossColor c.center.X (c.center.Y - 3) c.center.X (c.center.Y + 3) 1
+
+let drawCells (img: Image<Bgr, byte>) (drawCellContent: bool) (cells: Cell list) =
+    List.iter (fun c -> drawCell img drawCellContent c) cells
\ No newline at end of file
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Edges.fs b/Parasitemia/ParasitemiaCore/ImgTools/Edges.fs
new file mode 100644 (file)
index 0000000..b174ee9
--- /dev/null
@@ -0,0 +1,121 @@
+module ParasitemiaCore.Edges
+
+open System
+open System.Drawing
+open System.Collections.Generic
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+open Const
+open Histogram
+open Otsu
+
+/// <summary>
+/// Find edges of an image by using the Canny approach.
+/// The thresholds are automatically defined with otsu on gradient magnitudes.
+/// </summary>
+/// <param name="img"></param>
+let find (img: Image<Gray, float32>) : Matrix<byte> * Matrix<float32> * Matrix<float32> =
+    let w = img.Width
+    let h = img.Height
+
+    use sobelKernel =
+        new Matrix<float32>(array2D [[ 1.0f; 0.0f; -1.0f ]
+                                     [ 2.0f; 0.0f; -2.0f ]
+                                     [ 1.0f; 0.0f; -1.0f ]])
+
+    let xGradient = new Matrix<float32>(img.Size)
+    let yGradient = new Matrix<float32>(img.Size)
+    CvInvoke.Filter2D(img, xGradient, sobelKernel, Point(1, 1))
+    CvInvoke.Filter2D(img, yGradient, sobelKernel.Transpose(), Point(1, 1))
+
+    use magnitudes = new Matrix<float32>(xGradient.Size)
+    use angles = new Matrix<float32>(xGradient.Size)
+    CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles.
+
+    let thresholdHigh, thresholdLow =
+        let sensibilityHigh = 0.1f
+        let sensibilityLow = 0.0f
+        let threshold, _, _ = otsu (histogramMat magnitudes 300)
+        threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
+
+    // Non-maximum suppression.
+    use nms = new Matrix<byte>(xGradient.Size)
+
+    let nmsData = nms.Data
+    let anglesData = angles.Data
+    let magnitudesData = magnitudes.Data
+    let xGradientData = xGradient.Data
+    let yGradientData = yGradient.Data
+
+    for i in 0 .. h - 1 do
+        nmsData.[i, 0] <- 0uy
+        nmsData.[i, w - 1] <- 0uy
+
+    for j in 0 .. w - 1 do
+        nmsData.[0, j] <- 0uy
+        nmsData.[h - 1, j] <- 0uy
+
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            let vx = xGradientData.[i, j]
+            let vy = yGradientData.[i, j]
+            if vx <> 0.f || vy <> 0.f
+            then
+                let angle = anglesData.[i, j]
+
+                let vx', vy' = abs vx, abs vy
+                let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
+                let ratio1 = 1.f - ratio2
+
+                let mNeigbors (sign: int) : float32 =
+                    if angle < PI / 4.f
+                    then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
+                    elif angle < PI / 2.f
+                    then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j]
+                    elif angle < 3.f * PI / 4.f
+                    then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign]
+                    elif angle < PI
+                    then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign]
+                    elif angle < 5.f * PI / 4.f
+                    then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign]
+                    elif angle < 3.f * PI / 2.f
+                    then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j]
+                    elif angle < 7.f * PI / 4.f
+                    then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign]
+                    else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]
+
+                let m = magnitudesData.[i, j]
+                if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
+                then
+                    nmsData.[i, j] <- 1uy
+
+    // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
+
+    let edges = new Matrix<byte>(xGradient.Size)
+    let edgesData = edges.Data
+
+    // Hysteresis thresholding.
+    let toVisit = Stack<Point>()
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh
+            then
+                nmsData.[i, j] <- 0uy
+                toVisit.Push(Point(j, i))
+                while toVisit.Count > 0 do
+                    let p = toVisit.Pop()
+                    edgesData.[p.Y, p.X] <- 1uy
+                    for i' in -1 .. 1  do
+                        for j' in -1 .. 1 do
+                            if i' <> 0 || j' <> 0
+                            then
+                                let ni = p.Y + i'
+                                let nj = p.X + j'
+                                if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
+                                then
+                                    nmsData.[ni, nj] <- 0uy
+                                    toVisit.Push(Point(nj, ni))
+
+    edges, xGradient, yGradient
\ No newline at end of file
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Histogram.fs b/Parasitemia/ParasitemiaCore/ImgTools/Histogram.fs
new file mode 100644 (file)
index 0000000..92b73f6
--- /dev/null
@@ -0,0 +1,84 @@
+module ParasitemiaCore.Histogram
+
+open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+type Histogram = {
+    data: int[]
+    total: int // Number of elements.
+    sum: int // Sum of all intensity.
+    min: float32
+    max: float32 }
+
+let histogramImg (img: Image<Gray, float32>) (nbSamples: int) : Histogram =
+    let imgData = img.Data
+
+    let min, max =
+        let min = ref [| 0.0 |]
+        let minLocation = ref <| [| Point() |]
+        let max = ref [| 0.0 |]
+        let maxLocation = ref <| [| Point() |]
+        img.MinMax(min, max, minLocation, maxLocation)
+        float32 (!min).[0], float32 (!max).[0]
+
+    let inline bin (x: float32) : int =
+        let p = int ((x - min) / (max - min) * float32 nbSamples)
+        if p >= nbSamples then nbSamples - 1 else p
+
+    let data = Array.zeroCreate nbSamples
+
+    for i in 0 .. img.Height - 1 do
+        for j in 0 .. img.Width - 1 do
+            let p = bin imgData.[i, j, 0]
+            data.[p] <- data.[p] + 1
+
+    { data = data; total = img.Height * img.Width; sum = Array.sum data; min = min; max = max }
+
+let histogramMat (mat: Matrix<float32>) (nbSamples: int) : Histogram =
+    let matData = mat.Data
+
+    let min, max =
+        let min = ref 0.0
+        let minLocation = ref <| Point()
+        let max = ref 0.0
+        let maxLocation = ref <| Point()
+        mat.MinMax(min, max, minLocation, maxLocation)
+        float32 !min, float32 !max
+
+    let inline bin (x: float32) : int =
+        let p = int ((x - min) / (max - min) * float32 nbSamples)
+        if p >= nbSamples then nbSamples - 1 else p
+
+    let data = Array.zeroCreate nbSamples
+
+    for i in 0 .. mat.Height - 1 do
+        for j in 0 .. mat.Width - 1 do
+            let p = bin matData.[i, j]
+            data.[p] <- data.[p] + 1
+
+    { data = data; total = mat.Height * mat.Width; sum = Array.sum data; min = min; max = max }
+
+let histogram (values: float32 seq) (nbSamples: int) : Histogram =
+    let mutable min = Single.MaxValue
+    let mutable max = Single.MinValue
+    let mutable n = 0
+
+    for v in values do
+        n <- n + 1
+        if v < min then min <- v
+        if v > max then max <- v
+
+    let inline bin (x: float32) : int =
+        let p = int ((x - min) / (max - min) * float32 nbSamples)
+        if p >= nbSamples then nbSamples - 1 else p
+
+    let data = Array.zeroCreate nbSamples
+
+    for v in values do
+        let p = bin v
+        data.[p] <- data.[p] + 1
+
+    { data = data; total = n; sum = Array.sum data; min = min; max = max }
\ No newline at end of file
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/IO.fs b/Parasitemia/ParasitemiaCore/ImgTools/IO.fs
new file mode 100644 (file)
index 0000000..219eab6
--- /dev/null
@@ -0,0 +1,15 @@
+module ParasitemiaCore.IO
+
+open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+let saveImg (img: Image<'TColor, 'TDepth>) (filepath: string) =
+    img.Save(filepath)
+
+let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
+    use img = new Image<Gray, 'TDeph>(mat.Size)
+    mat.CopyTo(img)
+    saveImg img filepath
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/ImgTools.fs b/Parasitemia/ParasitemiaCore/ImgTools/ImgTools.fs
new file mode 100644 (file)
index 0000000..0510bc8
--- /dev/null
@@ -0,0 +1,49 @@
+module ParasitemiaCore.ImgTools
+
+open System
+open System.Drawing
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+let normalize (img: Image<Gray, float32>) (upperLimit: float) : Image<Gray, float32> =
+    let min = ref [| 0.0 |]
+    let minLocation = ref <| [| Point() |]
+    let max = ref [| 0.0 |]
+    let maxLocation = ref <| [| Point() |]
+    img.MinMax(min, max, minLocation, maxLocation)
+    let normalized = (img - (!min).[0]) / ((!max).[0] - (!min).[0])
+    if upperLimit = 1.0
+    then normalized
+    else upperLimit * normalized
+
+let mergeChannels (img: Image<Bgr, float32>) (rgbWeights: float * float * float) : Image<Gray, float32> =
+    match rgbWeights with
+    | 1., 0., 0. -> img.[2]
+    | 0., 1., 0. -> img.[1]
+    | 0., 0., 1. -> img.[0]
+    | redFactor, greenFactor, blueFactor ->
+        let result = new Image<Gray, float32>(img.Size)
+        CvInvoke.AddWeighted(result, 1., img.[2], redFactor, 0., result)
+        CvInvoke.AddWeighted(result, 1., img.[1], greenFactor, 0., result)
+        CvInvoke.AddWeighted(result, 1., img.[0], blueFactor, 0., result)
+        result
+
+let mergeChannelsWithProjection (img: Image<Bgr, float32>) (v1r: float32, v1g: float32, v1b: float32) (v2r: float32, v2g: float32, v2b: float32) (upperLimit: float)  : Image<Gray, float32> =
+    let vr, vg, vb = v2r - v1r, v2g - v1g, v2b - v1b
+    let vMagnitude = sqrt (vr ** 2.f + vg ** 2.f + vb ** 2.f)
+    let project (r: float32) (g: float32) (b: float32) = ((r - v1r) * vr + (g - v1g) * vg + (b - v1b) * vb) / vMagnitude
+    let result = new Image<Gray, float32>(img.Size)
+    // TODO: Essayer en bindant Data pour gagner du temps
+    for i in 0 .. img.Height - 1 do
+        for j in 0 .. img.Width - 1 do
+            result.Data.[i, j, 0] <- project img.Data.[i, j, 2] img.Data.[i, j, 1] img.Data.[i, j, 0]
+    normalize result upperLimit
+
+// Normalize image values between 0uy and 255uy.
+let normalizeAndConvert (img: Image<Gray, 'TDepth>) : Image<Gray, byte> =
+    (normalize (img.Convert<Gray, float32>()) 255.).Convert<Gray, byte>()
+
+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)
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs b/Parasitemia/ParasitemiaCore/ImgTools/Morpho.fs
new file mode 100644 (file)
index 0000000..b3a2e75
--- /dev/null
@@ -0,0 +1,670 @@
+module ParasitemiaCore.Morpho
+
+open System
+open System.Drawing
+open System.Collections.Generic
+open System.Linq
+
+open Emgu.CV
+open Emgu.CV.Structure
+
+open Types
+
+/// <summary>
+/// Remove M-adjacent pixels. It may be used after thinning.
+/// </summary>
+let suppressMAdjacency (img: Matrix<byte>) =
+    let w = img.Width
+    let h = img.Height
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            if img.[i, j] > 0uy && img.Data.[i + 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i - 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i - 1, j - 1] = 0uy)
+            then
+                img.[i, j] <- 0uy
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            if img.[i, j] > 0uy && img.Data.[i - 1, j] > 0uy && (img.Data.[i, j - 1] > 0uy && img.Data.[i + 1, j + 1] = 0uy || img.Data.[i, j + 1] > 0uy && img.Data.[i + 1, j - 1] = 0uy)
+            then
+                img.[i, j] <- 0uy
+
+type ExtremumType =
+    | Maxima = 1
+    | Minima = 2
+
+let findExtremum (img: Image<Gray, 'TDepth>) (extremumType: ExtremumType) : IEnumerable<Points> =
+    let w = img.Width
+    let h = img.Height
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
+
+    let imgData = img.Data
+    let suppress: bool[,] = Array2D.zeroCreate h w
+
+    let result = List<List<Point>>()
+
+    let flood (start: Point) : List<List<Point>> =
+        let sameLevelToCheck = Stack<Point>()
+        let betterLevelToCheck = Stack<Point>()
+        betterLevelToCheck.Push(start)
+
+        let result' = List<List<Point>>()
+
+        while betterLevelToCheck.Count > 0 do
+            let p = betterLevelToCheck.Pop()
+            if not suppress.[p.Y, p.X]
+            then
+                suppress.[p.Y, p.X] <- true
+                sameLevelToCheck.Push(p)
+                let current = List<Point>()
+
+                let mutable betterExists = false
+
+                while sameLevelToCheck.Count > 0 do
+                    let p' = sameLevelToCheck.Pop()
+                    let currentLevel = imgData.[p'.Y, p'.X, 0]
+                    current.Add(p') |> ignore
+                    for i, j in se do
+                        let ni = i + p'.Y
+                        let nj = j + p'.X
+                        if ni >= 0 && ni < h && nj >= 0 && nj < w
+                        then
+                            let level = imgData.[ni, nj, 0]
+                            let notSuppressed = not suppress.[ni, nj]
+
+                            if level = currentLevel && notSuppressed
+                            then
+                                suppress.[ni, nj] <- true
+                                sameLevelToCheck.Push(Point(nj, ni))
+                            elif if extremumType = ExtremumType.Maxima then level > currentLevel else level < currentLevel
+                            then
+                                betterExists <- true
+                                if notSuppressed
+                                then
+                                    betterLevelToCheck.Push(Point(nj, ni))
+
+                if not betterExists
+                then
+                    result'.Add(current)
+        result'
+
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            let maxima = flood (Point(j, i))
+            if maxima.Count > 0
+            then
+                result.AddRange(maxima)
+
+    result.Select(fun l -> Points(l))
+
+let findMaxima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
+    findExtremum img ExtremumType.Maxima
+
+let findMinima (img: Image<Gray, 'TDepth>) : IEnumerable<Points> =
+    findExtremum img ExtremumType.Minima
+
+type PriorityQueue () =
+    let size = 256
+    let q: Points[] = Array.init size (fun i -> Points())
+    let mutable highest = -1 // Value of the first elements of 'q'.
+    let mutable lowest = size
+
+    member this.NextMax () : byte * Point =
+        if this.IsEmpty
+        then
+            invalidOp "Queue is empty"
+        else
+            let l = q.[highest]
+            let next = l.First()
+            l.Remove(next) |> ignore
+            let value = byte highest
+
+            if l.Count = 0
+            then
+                highest <- highest - 1
+                while highest > lowest && q.[highest].Count = 0 do
+                    highest <- highest - 1
+                if highest = lowest
+                then
+                    highest <- -1
+                    lowest <- size
+
+            value, next
+
+    member this.NextMin () : byte * Point =
+        if this.IsEmpty
+        then
+            invalidOp "Queue is empty"
+        else
+            let l = q.[lowest + 1]
+            let next = l.First()
+            l.Remove(next) |> ignore
+            let value = byte (lowest + 1)
+
+            if l.Count = 0
+            then
+                lowest <- lowest + 1
+                while lowest < highest && q.[lowest + 1].Count = 0 do
+                    lowest <- lowest + 1
+                if highest = lowest
+                then
+                    highest <- -1
+                    lowest <- size
+
+            value, next
+
+    member this.Max =
+        highest |> byte
+
+    member this.Min =
+        lowest + 1 |> byte
+
+    member this.Add (value: byte) (p: Point) =
+        let vi = int value
+
+        if vi > highest
+        then
+            highest <- vi
+        if vi <= lowest
+        then
+            lowest <- vi - 1
+
+        q.[vi].Add(p) |> ignore
+
+    member this.Remove (value: byte) (p: Point) =
+        let vi = int value
+        if q.[vi].Remove(p) && q.[vi].Count = 0
+        then
+            if vi = highest
+            then
+                highest <- highest - 1
+                while highest > lowest && q.[highest].Count = 0 do
+                    highest <- highest - 1
+            elif vi - 1 = lowest
+            then
+                lowest <- lowest + 1
+                while lowest < highest && q.[lowest + 1].Count = 0 do
+                    lowest <- lowest + 1
+
+            if highest = lowest // The queue is now empty.
+            then
+                highest <- -1
+                lowest <- size
+
+    member this.IsEmpty =
+        highest = -1
+
+    member this.Clear () =
+        while highest > lowest  do
+            q.[highest].Clear()
+            highest <- highest - 1
+        highest <- -1
+        lowest <- size
+
+type private AreaState =
+    | Removed = 1
+    | Unprocessed = 2
+    | Validated = 3
+
+type private AreaOperation =
+    | Opening = 1
+    | Closing = 2
+
+[<AllowNullLiteral>]
+type private Area (elements: Points) =
+    member this.Elements = elements
+    member val Intensity = None with get, set
+    member val State = AreaState.Unprocessed with get, set
+
+let private areaOperation (img: Image<Gray, byte>) (area: int) (op: AreaOperation) =
+    let w = img.Width
+    let h = img.Height
+    let imgData = img.Data
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
+
+    let areas = List<Area>((if op = AreaOperation.Opening then findMaxima img else findMinima img) |> Seq.map Area)
+
+    let pixels: Area[,] = Array2D.create h w null
+    for m in areas do
+        for e in m.Elements do
+            pixels.[e.Y, e.X] <- m
+
+    let queue = PriorityQueue()
+
+    let addEdgeToQueue (elements: Points) =
+        for p in elements do
+            for i, j in se do
+                let ni = i + p.Y
+                let nj = j + p.X
+                let p' = Point(nj, ni)
+                if ni >= 0 && ni < h && nj >= 0 && nj < w && not (elements.Contains(p'))
+                then
+                    queue.Add (imgData.[ni, nj, 0]) p'
+
+    // Reverse order is quicker.
+    for i in areas.Count - 1 .. -1 .. 0 do
+        let m = areas.[i]
+        if m.Elements.Count <= area && m.State <> AreaState.Removed
+        then
+            queue.Clear()
+            addEdgeToQueue m.Elements
+
+            let mutable intensity = if op = AreaOperation.Opening then queue.Max else queue.Min
+            let nextElements = Points()
+
+            let mutable stop = false
+            while not stop do
+                let intensity', p = if op = AreaOperation.Opening then queue.NextMax () else queue.NextMin ()
+                let mutable merged = false
+
+                if intensity' = intensity // The intensity doesn't change.
+                then
+                    if m.Elements.Count + nextElements.Count + 1 > area
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity
+                        stop <- true
+                    else
+                        nextElements.Add(p) |> ignore
+
+                elif if op = AreaOperation.Opening then intensity' < intensity else intensity' > intensity
+                then
+                    m.Elements.UnionWith(nextElements)
+                    for e in nextElements do
+                        pixels.[e.Y, e.X] <- m
+
+                    if m.Elements.Count = area
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity')
+                        stop <- true
+                    else
+                        intensity <- intensity'
+                        nextElements.Clear()
+                        nextElements.Add(p) |> ignore
+
+                else
+                    match pixels.[p.Y, p.X] with
+                    | null -> ()
+                    | m' ->
+                        if m'.Elements.Count + m.Elements.Count <= area
+                        then
+                            m'.State <- AreaState.Removed
+                            for e in m'.Elements do
+                                pixels.[e.Y, e.X] <- m
+                                queue.Remove imgData.[e.Y, e.X, 0] e
+                            addEdgeToQueue m'.Elements
+                            m.Elements.UnionWith(m'.Elements)
+                            let intensityMax = if op = AreaOperation.Opening then queue.Max else queue.Min
+                            if intensityMax <> intensity
+                            then
+                                intensity <- intensityMax
+                                nextElements.Clear()
+                            merged <- true
+
+                    if not merged
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some (intensity)
+                        stop <- true
+
+                if not stop && not merged
+                then
+                    for i, j in se do
+                        let ni = i + p.Y
+                        let nj = j + p.X
+                        let p' = Point(nj, ni)
+                        if ni < 0 || ni >= h || nj < 0 || nj >= w
+                        then
+                            m.State <- AreaState.Validated
+                            m.Intensity <- Some (intensity)
+                            stop <- true
+                        elif not (m.Elements.Contains(p')) && not (nextElements.Contains(p'))
+                        then
+                            queue.Add (imgData.[ni, nj, 0]) p'
+
+                if queue.IsEmpty
+                then
+                    if m.Elements.Count + nextElements.Count <= area
+                    then
+                        m.State <- AreaState.Validated
+                        m.Intensity <- Some intensity'
+                        m.Elements.UnionWith(nextElements)
+                    stop <- true
+
+    for m in areas do
+        if m.State = AreaState.Validated
+        then
+            match m.Intensity with
+            | Some i ->
+                for p in m.Elements do
+                    imgData.[p.Y, p.X, 0] <- i
+            | _ -> ()
+    ()
+
+/// <summary>
+/// Area opening on byte image.
+/// </summary>
+let areaOpen (img: Image<Gray, byte>) (area: int) =
+    areaOperation img area AreaOperation.Opening
+
+/// <summary>
+/// Area closing on byte image.
+/// </summary>
+let areaClose (img: Image<Gray, byte>) (area: int) =
+    areaOperation img area AreaOperation.Closing
+
+// A simpler algorithm than 'areaOpen' on byte image but slower.
+let areaOpen2 (img: Image<Gray, byte>) (area: int) =
+    let w = img.Width
+    let h = img.Height
+    let imgData = img.Data
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
+
+    let histogram = Array.zeroCreate 256
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            let v = imgData.[i, j, 0] |> int
+            histogram.[v] <- histogram.[v] + 1
+
+    let flooded : bool[,] = Array2D.zeroCreate h w
+
+    let pointsChecked = HashSet<Point>()
+    let pointsToCheck = Stack<Point>()
+
+    for level in 255 .. -1 .. 0 do
+        let mutable n = histogram.[level]
+        if n > 0
+        then
+            for i in 0 .. h - 1 do
+                for j in 0 .. w - 1 do
+                    if not flooded.[i, j] && imgData.[i, j, 0] = byte level
+                    then
+                        let mutable maxNeighborValue = 0uy
+                        pointsChecked.Clear()
+                        pointsToCheck.Clear()
+                        pointsToCheck.Push(Point(j, i))
+
+                        while pointsToCheck.Count > 0 do
+                            let next = pointsToCheck.Pop()
+                            pointsChecked.Add(next) |> ignore
+                            flooded.[next.Y, next.X] <- true
+
+                            for nx, ny in se do
+                                let p = Point(next.X + nx, next.Y + ny)
+                                if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h
+                                then
+                                    let v = imgData.[p.Y, p.X, 0]
+                                    if v = byte level
+                                    then
+                                        if not (pointsChecked.Contains(p))
+                                        then
+                                            pointsToCheck.Push(p)
+                                    elif v > maxNeighborValue
+                                    then
+                                        maxNeighborValue <- v
+
+                        if int maxNeighborValue < level && pointsChecked.Count <= area
+                        then
+                            for p in pointsChecked do
+                                imgData.[p.Y, p.X, 0] <- maxNeighborValue
+
+[<AllowNullLiteral>]
+type Island (cmp: IComparer<float32>) =
+    member val Shore = Heap.Heap<float32, Point>(cmp) with get
+    member val Level = 0.f with get, set
+    member val Surface = 0 with get, set
+    member this.IsInfinite = this.Surface = Int32.MaxValue
+
+let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: ('a -> float32 -> unit) option) (op: AreaOperation) =
+    let w = img.Width
+    let h = img.Height
+    let earth = img.Data
+    let se = [| -1, 0; 0, -1; 1, 0; 0, 1 |]
+
+    let comparer = if op = AreaOperation.Opening
+                   then { new IComparer<float32> with member this.Compare(v1, v2) = v1.CompareTo(v2) }
+                   else { new IComparer<float32> with member this.Compare(v1, v2) = v2.CompareTo(v1) }
+
+    let ownership: Island[,] = Array2D.create h w null
+
+    // Initialize islands with their shore.
+    let islands = List<Island>()
+    let extremum = img |> if op = AreaOperation.Opening then findMaxima else findMinima
+    for e in extremum do
+        let island =
+            let p = e.First()
+            Island(comparer, Level = earth.[p.Y, p.X, 0], Surface = e.Count)
+        islands.Add(island)
+        let shorePoints = Points()
+        for p in e do
+            ownership.[p.Y, p.X] <- island
+            for i, j in se do
+                let ni = i + p.Y
+                let nj = j + p.X
+                let neighbor = Point(nj, ni)
+                if ni >= 0 && ni < h && nj >= 0 && nj < w && Object.ReferenceEquals(ownership.[ni, nj], null) && not (shorePoints.Contains(neighbor))
+                then
+                    shorePoints.Add(neighbor) |> ignore
+                    island.Shore.Add earth.[ni, nj, 0] neighbor
+
+    for area, obj in areas do
+        for island in islands do
+            let mutable stop = island.Shore.IsEmpty
+
+            // 'true' if 'p' is owned or adjacent to 'island'.
+            let inline ownedOrAdjacent (p: Point) : bool =
+                ownership.[p.Y, p.X] = island ||
+                (p.Y > 0 && ownership.[p.Y - 1, p.X] = island) ||
+                (p.Y < h - 1 && ownership.[p.Y + 1, p.X] = island) ||
+                (p.X > 0 && ownership.[p.Y, p.X - 1] = island) ||
+                (p.X < w - 1 && ownership.[p.Y, p.X + 1] = island)
+
+            while not stop && island.Surface < area do
+                let level, next = island.Shore.Max
+                let other = ownership.[next.Y, next.X]
+                if other = island // During merging, some points on the shore may be owned by the island itself -> ignored.
+                then
+                    island.Shore.RemoveNext ()
+                else
+                    if not <| Object.ReferenceEquals(other, null)
+                    then // We touching another island.
+                        if island.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area || comparer.Compare(island.Level, other.Level) < 0
+                        then
+                            stop <- true
+                        else // We can merge 'other' into 'surface'.
+                            island.Surface <- island.Surface + other.Surface
+                            island.Level <- other.Level
+                            // island.Level <- if comparer.Compare(island.Level, other.Level) > 0 then other.Level else island.Level
+                            for l, p in other.Shore do
+                                let mutable currentY = p.Y + 1
+                                while currentY < h && ownership.[currentY, p.X] = other do
+                                    ownership.[currentY, p.X] <- island
+                                    currentY <- currentY + 1
+                                island.Shore.Add l p
+                            other.Shore.Clear()
+
+                    elif comparer.Compare(level, island.Level) > 0
+                    then
+                        stop <- true
+                    else
+                        island.Shore.RemoveNext ()
+                        for i, j in se do
+                            let ni = i + next.Y
+                            let nj = j + next.X
+                            if ni < 0 || ni >= h || nj < 0 || nj >= w
+                            then
+                                island.Surface <- Int32.MaxValue
+                                stop <- true
+                            else
+                                let neighbor = Point(nj, ni)
+                                if not <| ownedOrAdjacent neighbor
+                                then
+                                    island.Shore.Add earth.[ni, nj, 0] neighbor
+                        if not stop
+                        then
+                            ownership.[next.Y, next.X] <- island
+                            island.Level <- level
+                            island.Surface <- island.Surface + 1
+
+        let mutable diff = 0.f
+
+        for i in 0 .. h - 1 do
+            for j in 0 .. w - 1 do
+                match ownership.[i, j] with
+                | null -> ()
+                | island ->
+                    let l = island.Level
+                    diff <- diff + l - earth.[i, j, 0]
+                    earth.[i, j, 0] <- l
+
+        match f with
+        | Some f' -> f' obj diff
+        | _ -> ()
+    ()
+
+/// <summary>
+/// Area opening on float image.
+/// </summary>
+let areaOpenF (img: Image<Gray, float32>) (area: int) =
+    areaOperationF img [ area, () ] None AreaOperation.Opening
+
+/// <summary>
+/// Area closing on float image.
+/// </summary>
+let areaCloseF (img: Image<Gray, float32>) (area: int) =
+    areaOperationF img [ area, () ] None AreaOperation.Closing
+
+/// <summary>
+/// Area closing on float image with different areas. Given areas must be sorted increasingly.
+/// For each area the function 'f' is called with the associated area value of type 'a and the volume difference
+/// Between the previous and the current closing.
+/// </summary>
+let areaOpenFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
+    areaOperationF img areas (Some f) AreaOperation.Opening
+
+/// <summary>
+/// Same as 'areaOpenFWithFun' for closing operation.
+/// </summary>
+let areaCloseFWithFun (img: Image<Gray, float32>) (areas: (int * 'a) list) (f: 'a -> float32 -> unit) =
+    areaOperationF img areas (Some f) AreaOperation.Closing
+
+/// <summary>
+/// Zhang and Suen thinning algorithm.
+/// Modify 'mat' in place.
+/// </summary>
+let thin (mat: Matrix<byte>) =
+    let w = mat.Width
+    let h = mat.Height
+    let mutable data1 = mat.Data
+    let mutable data2 = Array2D.copy data1
+
+    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 p2 = if i = 0 then 0uy else data1.[i-1, j]
+                    let p3 = if i = 0 || j = w-1 then 0uy else data1.[i-1, j+1]
+                    let p4 = if j = w-1 then 0uy else data1.[i, j+1]
+                    let p5 = if i = h-1 || j = w-1 then 0uy else data1.[i+1, j+1]
+                    let p6 = if i = h-1 then 0uy else data1.[i+1, j]
+                    let p7 = if i = h-1 || j = 0 then 0uy else data1.[i+1, j-1]
+                    let p8 = if j = 0 then 0uy else data1.[i, j-1]
+                    let p9 = if i = 0 || j = 0 then 0uy else data1.[i-1, j-1]
+
+                    let sumNeighbors = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
+                    if sumNeighbors >= 2uy && sumNeighbors <= 6uy &&
+                        (if p2 = 0uy && p3 = 1uy then 1 else 0) +
+                        (if p3 = 0uy && p4 = 1uy then 1 else 0) +
+                        (if p4 = 0uy && p5 = 1uy then 1 else 0) +
+                        (if p5 = 0uy && p6 = 1uy then 1 else 0) +
+                        (if p6 = 0uy && p7 = 1uy then 1 else 0) +
+                        (if p7 = 0uy && p8 = 1uy then 1 else 0) +
+                        (if p8 = 0uy && p9 = 1uy then 1 else 0) +
+                        (if p9 = 0uy && p2 = 1uy then 1 else 0) = 1 &&
+                        if oddIteration
+                        then p2 * p4 * p6 = 0uy && p4 * p6 * p8 = 0uy
+                        else p2 * p4 * p8 = 0uy && p2 * p6 * p8 = 0uy
+                    then
+                        data2.[i, j] <- 0uy
+                        pixelChanged <- true
+                else
+                    data2.[i, j] <- 0uy
+
+        oddIteration <- not oddIteration
+        let tmp = data1
+        data1 <- data2
+        data2 <- tmp
+
+/// <summary>
+/// Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
+/// Modify 'mat' in place.
+/// </summary>
+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
+
+    use 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<Point>()
+                let neighborsToCheck = Stack<Point>()
+                neighborsToCheck.Push(Point(j, i))
+                data'.[i, j] <- 0uy
+
+                while neighborsToCheck.Count > 0 do
+                    let n = neighborsToCheck.Pop()
+                    neighborhood.Add(n)
+                    for (ni, nj) in neighbors do
+                        let pi = n.Y + ni
+                        let pj = n.X + nj
+                        if pi >= 0 && pi < h && pj >= 0 && pj < w && data'.[pi, pj] = 1uy
+                        then
+                            neighborsToCheck.Push(Point(pj, pi))
+                            data'.[pi, pj] <- 0uy
+                if neighborhood.Count <= areaSize
+                then
+                    for n in neighborhood do
+                        data.[n.Y, n.X] <- 0uy
+
+let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Points =
+    let w = img.Width
+    let h = img.Height
+
+    let pointChecked = Points()
+    let pointToCheck = Stack<Point>(startPoints);
+
+    let data = img.Data
+
+    while pointToCheck.Count > 0 do
+        let next = pointToCheck.Pop()
+        pointChecked.Add(next) |> ignore
+        for ny in -1 .. 1 do
+            for nx in -1 .. 1 do
+                if ny <> 0 && nx <> 0
+                then
+                    let p = Point(next.X + nx, next.Y + ny)
+                    if p.X >= 0 && p.X < w && p.Y >= 0 && p.Y < h && data.[p.Y, p.X, 0] > 0uy && not (pointChecked.Contains p)
+                    then
+                        pointToCheck.Push(p)
+
+    pointChecked
diff --git a/Parasitemia/ParasitemiaCore/ImgTools/Otsu.fs b/Parasitemia/ParasitemiaCore/ImgTools/Otsu.fs
new file mode 100644 (file)
index 0000000..6b8ee53
--- /dev/null
@@ -0,0 +1,47 @@
+module ParasitemiaCore.Otsu
+
+open Histogram
+
+let otsu (hist: Histogram) : float32 * float32 * float32 =
+    let mutable sumB = 0
+    let mutable wB = 0
+    let mutable maximum = 0.0
+    let mutable level = 0
+    let sum = hist.data |> Array.mapi (fun i v -> i * v |> float) |> Array.sum
+
+    for i in 0 .. hist.data.Length - 1 do
+        wB <- wB + hist.data.[i]
+        if wB <> 0
+        then
+            let wF = hist.total - wB
+            if wF <> 0
+            then
+                sumB <- sumB + i * hist.data.[i]
+                let mB = (float sumB) / (float wB)
+                let mF = (sum - float sumB) / (float wF)
+                let between = (float wB) * (float wF) * (mB - mF) ** 2.;
+                if between >= maximum
+                then
+                    level <- i
+                    maximum <- between
+
+    let mean1 =
+        let mutable sum = 0
+        let mutable nb = 0
+        for i in 0 .. level - 1 do
+            sum <- sum + i * hist.data.[i]
+            nb <- nb + hist.data.[i]
+        (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2)
+
+    let mean2 =
+        let mutable sum = 0
+        let mutable nb = 0
+        for i in level + 1 .. hist.data.Length - 1 do
+            sum <- sum + i * hist.data.[i]
+            nb <- nb + hist.data.[i]
+        (sum + level * hist.data.[level] / 2) / (nb + hist.data.[level] / 2)
+
+    let toFloat l =
+        float32 l / float32 hist.data.Length * (hist.max - hist.min) + hist.min
+
+    toFloat level, toFloat mean1, toFloat mean2
\ No newline at end of file
index 9c85c7b..1d2ced3 100644 (file)
@@ -12,6 +12,7 @@ open Emgu.CV.Structure
 open Logger
 
 open Utils
+open Morpho
 open ImgTools
 open Config
 open Types
@@ -88,7 +89,7 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) (reportPr
             logTimeWithName "Parasites segmentation" (fun () -> reportWithVal 40 (ParasitesMarker.find img_parasites_filtered config))
 
         let! edges, xGradient, yGradient = logTimeWithName "Finding edges" (fun () ->
-            let edges, xGradient, yGradient = findEdges img_RBC_filtered
+            let edges, xGradient, yGradient = Edges.find img_RBC_filtered
             removeArea edges (config.RBCRadius.Pixel ** 2.f / 50.f |> int)
             reportWithVal 50 (edges, xGradient, yGradient))
 
@@ -109,43 +110,43 @@ let doAnalysis (img: Image<Bgr, byte>) (name: string) (config: Config) (reportPr
 
                 let buildFileName postfix = System.IO.Path.Combine(dirPath, name + postfix)
 
-                saveMat (edges * 255.0) (buildFileName " - edges.png")
+                IO.saveMat (edges * 255.0) (buildFileName " - edges.png")
 
-                saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png")
-                saveImg parasites.parasite (buildFileName " - parasites - stain.png")
-                saveImg parasites.nucleus (buildFileName " - parasites - infection.png")
+                IO.saveImg parasites.darkStain (buildFileName " - parasites - dark stain.png")
+                IO.saveImg parasites.parasite (buildFileName " - parasites - stain.png")
+                IO.saveImg parasites.nucleus (buildFileName " - parasites - infection.png")
 
                 let imgAllEllipses = img.Copy()
-                drawEllipses imgAllEllipses matchingEllipses.Ellipses (Bgr(255.0, 255.0, 255.0)) 0.04
-                saveImg imgAllEllipses (buildFileName " - ellipses - all.png")
+                Drawing.drawEllipses imgAllEllipses matchingEllipses.Ellipses (Bgr(255.0, 255.0, 255.0)) 0.04
+                IO.saveImg imgAllEllipses (buildFileName " - ellipses - all.png")
 
                 let imgEllipses = img_RBC_filtered.Convert<Bgr, byte>()
-                drawEllipses imgEllipses prunedEllipses (Bgr(0.0, 240.0, 240.0)) 1.0
-                saveImg imgEllipses (buildFileName " - ellipses.png")
+                Drawing.drawEllipses imgEllipses prunedEllipses (Bgr(0.0, 240.0, 240.0)) 1.0
+                IO.saveImg imgEllipses (buildFileName " - ellipses.png")
 
                 let imgCells = img.Copy()
-                drawCells imgCells false cells
-                saveImg imgCells (buildFileName " - cells.png")
+                Drawing.drawCells imgCells false cells
+                IO.saveImg imgCells (buildFileName " - cells.png")
 
                 let imgCells' = img.Copy()
-                drawCells imgCells' true cells
-                saveImg imgCells' (buildFileName " - cells - full.png")
+                Drawing.drawCells imgCells' true cells
+                IO.saveImg imgCells' (buildFileName " - cells - full.png")
 
                 let filteredGreenMaxima = gaussianFilter img_RBC config.LPFStandardDeviationRBC
                 for m in findMaxima filteredGreenMaxima do
-                    drawPoints filteredGreenMaxima m 255.f
-                saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png")
+                    Drawing.drawPoints filteredGreenMaxima m 255.f
+                IO.saveImg filteredGreenMaxima (buildFileName " - filtered - maxima.png")
 
-                saveImg img_RBC_filtered (buildFileName " - filtered.png")
-                saveImg imgWhitoutParasite (buildFileName " - filtered closed stain.png")
-                saveImg imgWithoutNucleus (buildFileName " - filtered closed infection.png")
+                IO.saveImg img_RBC_filtered (buildFileName " - filtered.png")
+                IO.saveImg imgWhitoutParasite (buildFileName " - filtered closed stain.png")
+                IO.saveImg imgWithoutNucleus (buildFileName " - filtered closed infection.png")
 
-                saveImg img_RBC (buildFileName " - source - RBC.png")
-                saveImg img_parasites (buildFileName " - source - parasites.png")
+                IO.saveImg img_RBC (buildFileName " - source - RBC.png")
+                IO.saveImg img_parasites (buildFileName " - source - parasites.png")
 
-                saveImg (normalize img_float.[2] 255.) (buildFileName " - source - red.png")
-                saveImg (normalize img_float.[1] 255.) (buildFileName " - source - green.png")
-                saveImg (normalize img_float.[0] 255.) (buildFileName " - source - blue.png")
+                IO.saveImg (normalize img_float.[2] 255.) (buildFileName " - source - red.png")
+                IO.saveImg (normalize img_float.[1] 255.) (buildFileName " - source - green.png")
+                IO.saveImg (normalize img_float.[0] 255.) (buildFileName " - source - blue.png")
             | _ -> ()
 
         return cells }
index 58886c0..6a9eb75 100644 (file)
     <Compile Include="Types.fs" />
     <Compile Include="EEOver.fs" />
     <Compile Include="Utils.fs" />
-    <Compile Include="ImgTools.fs" />
+    <Compile Include="ImgTools\Histogram.fs" />
+    <Compile Include="ImgTools\Otsu.fs" />
+    <Compile Include="ImgTools\Drawing.fs" />
+    <Compile Include="ImgTools\Edges.fs" />
+    <Compile Include="ImgTools\Morpho.fs" />
+    <Compile Include="ImgTools\IO.fs" />
+    <Compile Include="ImgTools\ImgTools.fs" />
     <Compile Include="Granulometry.fs" />
     <Compile Include="Config.fs" />
     <Compile Include="KMedians.fs" />
index 6c25d04..0d7c7ae 100644 (file)
@@ -7,6 +7,9 @@ open Emgu.CV
 open Emgu.CV.Structure
 
 open Utils
+open Histogram
+open Otsu
+open Morpho
 open ImgTools
 
 type Result = {
index 0d84d16..36c9495 100644 (file)
@@ -1,9 +1,9 @@
 module ParasitemiaCore.UnitsOfMeasure
 
-[<Measure>] type px
+[<Measure>] type px // Pixel.
 [<Measure>] type μm
 [<Measure>] type inch
-[<Measure>] type ppi = px / inch
+[<Measure>] type ppi = px / inch // Pixel per inch.
 
 let μmInchRatio = 25.4e3<μm/inch>
 
index 3cb0de4..e35d234 100644 (file)
@@ -3,7 +3,7 @@
         xmlns:d="http://schemas.microsoft.com/expression/blend/2008"
         xmlns:mc="http://schemas.openxmlformats.org/markup-compatibility/2006"
         mc:Ignorable="d"
-        x:Name="AboutWindow" Height="200.969" Width="282.313" MinHeight="100" MinWidth="100" Title="About" Icon="pack://application:,,,/Resources/icon.ico">
+        x:Name="AboutWindow" Height="220" Width="280" MinHeight="220" MinWidth="280" Title="About" Icon="pack://application:,,,/Resources/icon.ico" ResizeMode="NoResize">
    <Grid>
       <Grid.RowDefinitions>
          <RowDefinition Height="Auto"/>