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