Project the colors to have the best contrast for RBCs and parasites analyze.
[master-thesis.git] / Parasitemia / ParasitemiaCore / ImgTools.fs
index cc15af7..51b0c19 100644 (file)
@@ -13,6 +13,40 @@ 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> =
     let min = ref [| 0.0 |]
@@ -43,7 +77,7 @@ let histogramImg (img: Image<Gray, float32>) (nbSamples: int) : Histogram =
         img.MinMax(min, max, minLocation, maxLocation)
         float32 (!min).[0], float32 (!max).[0]
 
-    let bin (x: float32) : int =
+    let inline bin (x: float32) : int =
         let p = int ((x - min) / (max - min) * float32 nbSamples)
         if p >= nbSamples then nbSamples - 1 else p
 
@@ -67,7 +101,7 @@ let histogramMat (mat: Matrix<float32>) (nbSamples: int) : Histogram =
         mat.MinMax(min, max, minLocation, maxLocation)
         float32 !min, float32 !max
 
-    let bin (x: float32) : int =
+    let inline bin (x: float32) : int =
         let p = int ((x - min) / (max - min) * float32 nbSamples)
         if p >= nbSamples then nbSamples - 1 else p
 
@@ -90,7 +124,7 @@ let histogram (values: float32 seq) (nbSamples: int) : Histogram =
         if v < min then min <- v
         if v > max then max <- v
 
-    let bin (x: float32) : int =
+    let inline bin (x: float32) : int =
         let p = int ((x - min) / (max - min) * float32 nbSamples)
         if p >= nbSamples then nbSamples - 1 else p
 
@@ -107,7 +141,7 @@ let otsu (hist: Histogram) : float32 * float32 * float32 =
     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
+    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]
@@ -168,31 +202,19 @@ let suppressMAdjacency (img: Matrix<byte>) =
 /// The thresholds are automatically defined with otsu on gradient magnitudes.
 /// </summary>
 /// <param name="img"></param>
-let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32> * Image<Gray, float32> =
+let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Matrix<float32> * Matrix<float32> =
     let w = img.Width
     let h = img.Height
 
     use sobelKernel =
-        new ConvolutionKernelF(array2D [[ 1.0f; 0.0f; -1.0f ]
-                                        [ 2.0f; 0.0f; -2.0f ]
-                                        [ 1.0f; 0.0f; -1.0f ]], Point(1, 1))
-
-    let xGradient = img.Convolution(sobelKernel)
-    let yGradient = img.Convolution(sobelKernel.Transpose())
+        new Matrix<float32>(array2D [[ 1.0f; 0.0f; -1.0f ]
+                                     [ 2.0f; 0.0f; -2.0f ]
+                                     [ 1.0f; 0.0f; -1.0f ]])
 
-    let xGradientData = xGradient.Data
-    let yGradientData = yGradient.Data
-    for r in 0 .. h - 1 do
-        xGradientData.[r, 0, 0] <- 0.f
-        xGradientData.[r, w - 1, 0] <- 0.f
-        yGradientData.[r, 0, 0] <- 0.f
-        yGradientData.[r, w - 1, 0] <- 0.f
-
-    for c in 0 .. w - 1 do
-        xGradientData.[0, c, 0] <- 0.f
-        xGradientData.[h - 1, c, 0] <- 0.f
-        yGradientData.[0, c, 0] <- 0.f
-        yGradientData.[h - 1, c, 0] <- 0.f
+    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)
@@ -201,7 +223,6 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
     let thresholdHigh, thresholdLow =
         let sensibilityHigh = 0.1f
         let sensibilityLow = 0.0f
-        use magnitudesByte = magnitudes.Convert<byte>()
         let threshold, _, _ = otsu (histogramMat magnitudes 300)
         threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
 
@@ -224,8 +245,8 @@ let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float32>
 
     for i in 1 .. h - 2 do
         for j in 1 .. w - 2 do
-            let vx = xGradientData.[i, j, 0]
-            let vy = yGradientData.[i, j, 0]
+            let vx = xGradientData.[i, j]
+            let vy = yGradientData.[i, j]
             if vx <> 0.f || vy <> 0.f
             then
                 let angle = anglesData.[i, j]
@@ -733,12 +754,13 @@ let private areaOperationF (img: Image<Gray, float32>) (areas: (int * 'a) list)
                 else
                     if not <| Object.ReferenceEquals(other, null)
                     then // We touching another island.
-                        if island.IsInfinite || other.IsInfinite || island.Surface + other.Surface >= area
+                        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 <- if comparer.Compare(island.Level, other.Level) > 0 then island.Level else other.Level
+                            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
@@ -943,7 +965,7 @@ let drawLineF (img: Image<'TColor, 'TDepth>) (color: 'TColor) (x0: float) (y0: f
 let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Ellipse) (color: 'TColor) (alpha: float) =
     if alpha >= 1.0
     then
-        img.Draw(Emgu.CV.Structure.Ellipse(PointF(float32 e.Cx, float32 e.Cy), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
+        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))
@@ -957,7 +979,7 @@ let drawEllipse (img: Image<'TColor, 'TDepth>) (e: Ellipse) (color: 'TColor) (al
         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(float32 <| (e.A + 5.f + gapX) , float32 <| (e.A + 5.f + gapY)), SizeF(2.f * e.B, 2.f * e.A), e.Alpha / PI * 180.f), color, 1, CvEnum.LineType.AntiAlias)
+            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