1
module ParasitemiaCore.Edges
4 open System.Collections.Generic
13 // Sensibilities of the hysteresis search.
14 let sensibilityHigh = 0.1f
15 let sensibilityLow = 0.0f
18 /// Find edges of an image by using the Canny approach.
19 /// The thresholds are automatically defined with otsu on gradient magnitudes.
21 /// <param name="img"></param>
22 let find (img
: Image<Gray, float32
>) : Matrix<byte
> * Matrix<float32
> * Matrix<float32
> =
28 array2D
[[ -1.0f; 0.0f; 1.0f ]
30 [ -1.0f; 0.0f; 1.0f ]]
33 let xGradient = new Matrix<float32
> (img
.Size)
34 let yGradient = new Matrix<float32
> (img
.Size)
35 CvInvoke.Filter2D (img
, xGradient, sobelKernel, Point (1, 1))
36 CvInvoke.Filter2D (img
, yGradient, sobelKernel.Transpose (), Point (1, 1))
38 use magnitudes = new Matrix<float32
> (xGradient.Size)
39 use angles = new Matrix<float32
> (xGradient.Size)
40 CvInvoke.CartToPolar (xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles. The angles are between 0 and 2 * pi.
42 let thresholdHigh, thresholdLow
=
43 let threshold, _
, _
= otsu
(histogramMat
magnitudes 300)
44 threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
46 // Non-maximum suppression.
47 use nms = new Matrix<byte
> (xGradient.Size)
49 let nmsData = nms.Data
50 let anglesData = angles.Data
51 let magnitudesData = magnitudes.Data
52 let xGradientData = xGradient.Data
53 let yGradientData = yGradient.Data
57 let vx = xGradientData.[i
, j
]
58 let vy = yGradientData.[i
, j
]
59 if vx <> 0.f
|| vy <> 0.f
then
60 let angle = anglesData.[i
, j
]
62 let vx', vy' = abs
vx, abs
vy
63 let ratio2 = if vx' > vy' then vy' / vx' else vx' / vy'
64 let ratio1 = 1.f
- ratio2
66 let mNeigbors (sign
: int) : float32
=
67 if angle < PI / 4.f
then
68 ratio1 * magnitudesData.[i
, j
+ sign
] + ratio2 * magnitudesData.[i
+ sign
, j
+ sign
]
69 elif
angle < PI / 2.f
then
70 ratio2 * magnitudesData.[i
+ sign
, j
+ sign
] + ratio1 * magnitudesData.[i
+ sign
, j
]
71 elif
angle < 3.f
* PI / 4.f
then
72 ratio1 * magnitudesData.[i
+ sign
, j
] + ratio2 * magnitudesData.[i
+ sign
, j
- sign
]
74 ratio2 * magnitudesData.[i
+ sign
, j
- sign
] + ratio1 * magnitudesData.[i
, j
- sign
]
75 elif
angle < 5.f
* PI / 4.f
then
76 ratio1 * magnitudesData.[i
, j
- sign
] + ratio2 * magnitudesData.[i
- sign
, j
- sign
]
77 elif
angle < 3.f
* PI / 2.f
then
78 ratio2 * magnitudesData.[i
- sign
, j
- sign
] + ratio1 * magnitudesData.[i
- sign
, j
]
79 elif
angle < 7.f
* PI / 4.f
then
80 ratio1 * magnitudesData.[i
- sign
, j
] + ratio2 * magnitudesData.[i
- sign
, j
+ sign
]
82 ratio2 * magnitudesData.[i
- sign
, j
+ sign
] + ratio1 * magnitudesData.[i
, j
+ sign
]
84 let m = magnitudesData.[i
, j
]
85 if m >= thresholdLow
&& m > mNeigbors 1 && m > mNeigbors -1 then
88 // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
90 let edges = new Matrix<byte
> (xGradient.Size)
91 let edgesData = edges.Data
93 // Hysteresis thresholding.
94 let toVisit = Stack<Point> ()
97 if nmsData.[i
, j
] = 1uy && magnitudesData.[i
, j
] >= thresholdHigh then
99 toVisit.Push (Point (j
, i
))
100 while toVisit.Count > 0 do
101 let p = toVisit.Pop ()
102 edgesData.[p.Y, p.X] <- 1uy
105 if i
' <> 0 || j' <> 0 then
108 if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy then
109 nmsData.[ni, nj] <- 0uy
110 toVisit.Push (Point (nj, ni))
112 edges, xGradient, yGradient