Use k-means instead of k-medians.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index 2f021c8..5e32ccb 100644 (file)
@@ -21,6 +21,130 @@ let normalizeAndConvert (img: Image<Gray, float32>) : Image<Gray, byte> =
     ((img - (!min).[0]) / ((!max).[0] - (!min).[0]) * 255.0).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
+
+
+let suppressMConnections (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
+
+let findEdges (img: Image<Gray, float32>) : Matrix<byte> * Image<Gray, float> * Image<Gray, float> =
+    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).Convert<Gray, float>()
+    let yGradient = img.Convolution(sobelKernel.Transpose()).Convert<Gray, float>()
+
+    let xGradientData = xGradient.Data
+    let yGradientData = yGradient.Data
+    for r in 0 .. h - 1 do
+        xGradientData.[r, 0, 0] <- 0.0
+        xGradientData.[r, w - 1, 0] <- 0.0
+        yGradientData.[r, 0, 0] <- 0.0
+        yGradientData.[r, w - 1, 0] <- 0.0
+
+    for c in 0 .. w - 1 do
+        xGradientData.[0, c, 0] <- 0.0
+        xGradientData.[h - 1, c, 0] <- 0.0
+        yGradientData.[0, c, 0] <- 0.0
+        yGradientData.[h - 1, c, 0] <- 0.0
+
+    use magnitudes = new Matrix<float>(xGradient.Size)
+    CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, new Mat()) // Compute the magnitudes (without angles).
+
+    let thresholdHigh, thresholdLow =
+        let sensibility = 0.1
+        use magnitudesByte = magnitudes.Convert<byte>()
+        let threshold = CvInvoke.Threshold(magnitudesByte, magnitudesByte, 0.0, 1.0, CvEnum.ThresholdType.Otsu ||| CvEnum.ThresholdType.Binary)
+        threshold + (sensibility * threshold), threshold - (sensibility * threshold)
+
+    // Non-maximum suppression.
+    use nms = new Matrix<byte>(xGradient.Size)
+    nms.SetValue(1.0)
+
+    for i in 0 .. h - 1 do
+        nms.Data.[i, 0] <- 0uy
+        nms.Data.[i, w - 1] <- 0uy
+
+    for j in 0 .. w - 1 do
+        nms.Data.[0, j] <- 0uy
+        nms.Data.[h - 1, j] <- 0uy
+
+    for i in 1 .. h - 2 do
+        for j in 1 .. w - 2 do
+            let vx = xGradient.Data.[i, j, 0]
+            let vy = yGradient.Data.[i, j, 0]
+            let angle =
+                let a = atan2 vy vx
+                if a < 0.0 then 2. * Math.PI + a else a
+
+            let mNeigbors (sign: int) : float =
+                if angle < Math.PI / 8. || angle >= 15.0 * Math.PI / 8. then magnitudes.Data.[i, j + sign]
+                elif angle < 3.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j + sign]
+                elif angle < 5.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j]
+                elif angle < 7.0 * Math.PI / 8. then magnitudes.Data.[i + sign, j - sign]
+                elif angle < 9.0 * Math.PI / 8. then magnitudes.Data.[i, j - sign]
+                elif angle < 11.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j - sign]
+                elif angle < 13.0 * Math.PI / 8. then magnitudes.Data.[i - sign, j]
+                else magnitudes.Data.[i - sign, j + sign]
+
+            let m = magnitudes.Data.[i, j]
+            if m < mNeigbors 1 || m < mNeigbors -1 || m < thresholdLow
+            then
+                nms.Data.[i, j] <- 0uy
+
+    // suppressMConnections nms // It's not usefull for the rest of the process (ellipse detection).
+
+    let edges = new Matrix<byte>(xGradient.Size)
+
+    // Histeresis thresholding.
+    let toVisit = Stack<Point>()
+    for i in 0 .. h - 1 do
+        for j in 0 .. w - 1 do
+            if nms.Data.[i, j] = 1uy && magnitudes.Data.[i, j] >= thresholdHigh
+            then
+                nms.Data.[i, j] <- 0uy
+                toVisit.Push(Point(j, i))
+                while toVisit.Count > 0 do
+                    let p = toVisit.Pop()
+                    edges.Data.[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 && nms.Data.[ni, nj] = 1uy
+                                then
+                                    nms.Data.[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)
@@ -32,7 +156,6 @@ let drawPoints (img: Image<Gray, byte>) (points: Points) (intensity: byte) =
     for p in points do
         img.Data.[p.Y, p.X, 0] <- intensity
 
-
 type ExtremumType =
     | Maxima = 1
     | Minima = 2
@@ -356,6 +479,7 @@ let areaOpen (img: Image<Gray, byte>) (area: int) =
 let areaClose (img: Image<Gray, byte>) (area: int) =
     areaOperation img area AreaOperation.Closing
 
+// A simpler algorithm than 'areaOpen' but slower.
 let areaOpen2 (img: Image<Gray, byte>) (area: int) =
     let w = img.Width
     let h = img.Height
@@ -462,12 +586,6 @@ let thin (mat: Matrix<byte>) =
         data2 <- tmp
 
 
-// FIXME: replace by a queue or stack.
-let pop (l: List<'a>) : 'a =
-    let n = l.[l.Count - 1]
-    l.RemoveAt(l.Count - 1)
-    n
-
 // Remove all 8-connected pixels with an area equal or greater than 'areaSize'.
 // Modify 'mat' in place.
 let removeArea (mat: Matrix<byte>) (areaSize: int) =
@@ -493,37 +611,37 @@ let removeArea (mat: Matrix<byte>) (areaSize: int) =
         for j in 0..w-1 do
             if data'.[i, j] = 1uy
             then
-                let neighborhood = List<(int*int)>()
-                let neighborsToCheck = List<(int*int)>()
-                neighborsToCheck.Add((i, j))
+                let neighborhood = List<Point>()
+                let neighborsToCheck = Stack<Point>()
+                neighborsToCheck.Push(Point(j, i))
                 data'.[i, j] <- 0uy
 
                 while neighborsToCheck.Count > 0 do
-                    let (ci, cj) = pop neighborsToCheck
-                    neighborhood.Add((ci, cj))
+                    let n = neighborsToCheck.Pop()
+                    neighborhood.Add(n)
                     for (ni, nj) in neighbors do
-                        let pi = ci + ni
-                        let pj = cj + nj
+                        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.Add((pi, pj))
+                            neighborsToCheck.Push(Point(pj, pi))
                             data'.[pi, pj] <- 0uy
                 if neighborhood.Count <= areaSize
                 then
-                    for (ni, nj) in neighborhood do
-                        data.[ni, nj] <- 0uy
+                    for n in neighborhood do
+                        data.[n.Y, n.X] <- 0uy
 
 let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : List<Point> =
     let w = img.Width
     let h = img.Height
 
     let pointChecked = Points()
-    let pointToCheck = List<Point>(startPoints);
+    let pointToCheck = Stack<Point>(startPoints);
 
     let data = img.Data
 
     while pointToCheck.Count > 0 do
-        let next = pop pointToCheck
+        let next = pointToCheck.Pop()
         pointChecked.Add(next) |> ignore
         for ny in -1 .. 1 do
             for nx in -1 .. 1 do
@@ -532,26 +650,19 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
                     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.Add(p)
+                        pointToCheck.Push(p)
 
     List<Point>(pointChecked)
 
 
-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
-
 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: Types.Ellipse) (color: 'TColor) (alpha: float) =
 
     if alpha >= 1.0