* Treat some special cases when ellipses intersecting.
[master-thesis.git] / Parasitemia / Parasitemia / ImgTools.fs
index 2f021c8..f567e68 100644 (file)
@@ -21,6 +21,124 @@ 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 +150,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 +473,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
@@ -537,21 +655,14 @@ let connectedComponents (img: Image<Gray, byte>) (startPoints: List<Point>) : Li
     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