Split the module 'ImgTools' in many modules.
[master-thesis.git] / Parasitemia / ParasitemiaCore / ImgTools / Edges.fs
1 module ParasitemiaCore.Edges
2
3 open System
4 open System.Drawing
5 open System.Collections.Generic
6
7 open Emgu.CV
8 open Emgu.CV.Structure
9
10 open Const
11 open Histogram
12 open Otsu
13
14 /// <summary>
15 /// Find edges of an image by using the Canny approach.
16 /// The thresholds are automatically defined with otsu on gradient magnitudes.
17 /// </summary>
18 /// <param name="img"></param>
19 let find (img: Image<Gray, float32>) : Matrix<byte> * Matrix<float32> * Matrix<float32> =
20 let w = img.Width
21 let h = img.Height
22
23 use sobelKernel =
24 new Matrix<float32>(array2D [[ 1.0f; 0.0f; -1.0f ]
25 [ 2.0f; 0.0f; -2.0f ]
26 [ 1.0f; 0.0f; -1.0f ]])
27
28 let xGradient = new Matrix<float32>(img.Size)
29 let yGradient = new Matrix<float32>(img.Size)
30 CvInvoke.Filter2D(img, xGradient, sobelKernel, Point(1, 1))
31 CvInvoke.Filter2D(img, yGradient, sobelKernel.Transpose(), Point(1, 1))
32
33 use magnitudes = new Matrix<float32>(xGradient.Size)
34 use angles = new Matrix<float32>(xGradient.Size)
35 CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles.
36
37 let thresholdHigh, thresholdLow =
38 let sensibilityHigh = 0.1f
39 let sensibilityLow = 0.0f
40 let threshold, _, _ = otsu (histogramMat magnitudes 300)
41 threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
42
43 // Non-maximum suppression.
44 use nms = new Matrix<byte>(xGradient.Size)
45
46 let nmsData = nms.Data
47 let anglesData = angles.Data
48 let magnitudesData = magnitudes.Data
49 let xGradientData = xGradient.Data
50 let yGradientData = yGradient.Data
51
52 for i in 0 .. h - 1 do
53 nmsData.[i, 0] <- 0uy
54 nmsData.[i, w - 1] <- 0uy
55
56 for j in 0 .. w - 1 do
57 nmsData.[0, j] <- 0uy
58 nmsData.[h - 1, j] <- 0uy
59
60 for i in 1 .. h - 2 do
61 for j in 1 .. w - 2 do
62 let vx = xGradientData.[i, j]
63 let vy = yGradientData.[i, j]
64 if vx <> 0.f || vy <> 0.f
65 then
66 let angle = anglesData.[i, j]
67
68 let vx', vy' = abs vx, abs vy
69 let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
70 let ratio1 = 1.f - ratio2
71
72 let mNeigbors (sign: int) : float32 =
73 if angle < PI / 4.f
74 then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
75 elif angle < PI / 2.f
76 then ratio2 * magnitudesData.[i + sign, j + sign] + ratio1 * magnitudesData.[i + sign, j]
77 elif angle < 3.f * PI / 4.f
78 then ratio1 * magnitudesData.[i + sign, j] + ratio2 * magnitudesData.[i + sign, j - sign]
79 elif angle < PI
80 then ratio2 * magnitudesData.[i + sign, j - sign] + ratio1 * magnitudesData.[i, j - sign]
81 elif angle < 5.f * PI / 4.f
82 then ratio1 * magnitudesData.[i, j - sign] + ratio2 * magnitudesData.[i - sign, j - sign]
83 elif angle < 3.f * PI / 2.f
84 then ratio2 * magnitudesData.[i - sign, j - sign] + ratio1 * magnitudesData.[i - sign, j]
85 elif angle < 7.f * PI / 4.f
86 then ratio1 * magnitudesData.[i - sign, j] + ratio2 * magnitudesData.[i - sign, j + sign]
87 else ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]
88
89 let m = magnitudesData.[i, j]
90 if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
91 then
92 nmsData.[i, j] <- 1uy
93
94 // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
95
96 let edges = new Matrix<byte>(xGradient.Size)
97 let edgesData = edges.Data
98
99 // Hysteresis thresholding.
100 let toVisit = Stack<Point>()
101 for i in 0 .. h - 1 do
102 for j in 0 .. w - 1 do
103 if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh
104 then
105 nmsData.[i, j] <- 0uy
106 toVisit.Push(Point(j, i))
107 while toVisit.Count > 0 do
108 let p = toVisit.Pop()
109 edgesData.[p.Y, p.X] <- 1uy
110 for i' in -1 .. 1 do
111 for j' in -1 .. 1 do
112 if i' <> 0 || j' <> 0
113 then
114 let ni = p.Y + i'
115 let nj = p.X + j'
116 if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
117 then
118 nmsData.[ni, nj] <- 0uy
119 toVisit.Push(Point(nj, ni))
120
121 edges, xGradient, yGradient