* Try another approach to remove false ellipses without success (commented).

[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs

diff --git a/Parasitemia/Parasitemia/ImgTools.fs b/Parasitemia/Parasitemia/ImgTools.fs
index f64b1c3..3cfdc89 100644 (file)
--- a/Parasitemia/Parasitemia/ImgTools.fs
+++ b/Parasitemia/Parasitemia/ImgTools.fs
@@ -32,6 +32,123 @@ let saveMat (mat: Matrix<'TDepth>) (filepath: string) =
     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 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 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 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) |> Array.sum |> float
+
+    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
+
+
 let suppressMConnections (img: Matrix<byte>) =
     let w = img.Width
     let h = img.Height
@@ -79,9 +196,11 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
 
     let thresholdHigh, thresholdLow =
         let sensibilityHigh = 0.1f
-        let sensibilityLow = 0.1f
+        let sensibilityLow = 0.0f
         use magnitudesByte = magnitudes.Convert<byte>()
         let threshold = float32 <| CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
+        let threshold, _, _ = otsu (histogramMat magnitudes 300)
+
         threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
 
     // Non-maximum suppression.