1
module ParasitemiaCore.Edges
5 open System.Collections.Generic
14 // Sensibilities of the hysteresis search.
15 let sensibilityHigh = 0.1f
16 let sensibilityLow = 0.0f
19 /// Find edges of an image by using the Canny approach.
20 /// The thresholds are automatically defined with otsu on gradient magnitudes.
22 /// <param name="img"></param>
23 let find (img
: Image<Gray, float32
>) : Matrix<byte
> * Matrix<float32
> * Matrix<float32
> =
29 array2D
[[ -1.0f; 0.0f; 1.0f ]
31 [ -1.0f; 0.0f; 1.0f ]]
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))
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.
43 let thresholdHigh, thresholdLow
=
44 let threshold, _
, _
= otsu
(histogramMat
magnitudes 300)
45 threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
47 // Non-maximum suppression.
48 use nms = new Matrix<byte
> (xGradient.Size)
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
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
]
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
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
]
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
]
83 ratio2 * magnitudesData.[i
- sign
, j
+ sign
] + ratio1 * magnitudesData.[i
, j
+ sign
]
85 let m = magnitudesData.[i
, j
]
86 if m >= thresholdLow
&& m > mNeigbors 1 && m > mNeigbors -1 then
89 // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
91 let edges = new Matrix<byte
> (xGradient.Size)
92 let edgesData = edges.Data
94 // Hysteresis thresholding.
95 let toVisit = Stack<Point> ()
98 if nmsData.[i
, j
] = 1uy && magnitudesData.[i
, j
] >= thresholdHigh then
100 toVisit.Push (Point (j
, i
))
101 while toVisit.Count > 0 do
102 let p = toVisit.Pop ()
103 edgesData.[p.Y, p.X] <- 1uy
106 if i
' <> 0 || j' <> 0 then
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))
113 edges, xGradient, yGradient