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
> =
28 new Matrix<float32
>(array2D
[[ -1.0f; 0.0f; 1.0f ]
30 [ -1.0f; 0.0f; 1.0f ]])
32 let xGradient = new Matrix<float32
>(img
.Size)
33 let yGradient = new Matrix<float32
>(img
.Size)
34 CvInvoke.Filter2D(img
, xGradient, sobelKernel, Point(1, 1))
35 CvInvoke.Filter2D(img
, yGradient, sobelKernel.Transpose(), Point(1, 1))
37 use magnitudes = new Matrix<float32
>(xGradient.Size)
38 use angles = new Matrix<float32
>(xGradient.Size)
39 CvInvoke.CartToPolar(xGradient, yGradient, magnitudes, angles) // Compute the magnitudes and angles. The angles are between 0 and 2 * pi.
41 let thresholdHigh, thresholdLow
=
42 let threshold, _
, _
= otsu
(histogramMat
magnitudes 300)
43 threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
45 // Non-maximum suppression.
46 use nms = new Matrix<byte
>(xGradient.Size)
48 let nmsData = nms.Data
49 let anglesData = angles.Data
50 let magnitudesData = magnitudes.Data
51 let xGradientData = xGradient.Data
52 let yGradientData = yGradient.Data
54 for i
in 1 .. h - 2 do
55 for j
in 1 .. w - 2 do
56 let vx = xGradientData.[i
, j
]
57 let vy = yGradientData.[i
, j
]
58 if vx <> 0.f
|| vy <> 0.f
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
=
68 then ratio1 * magnitudesData.[i
, j
+ sign
] + ratio2 * magnitudesData.[i
+ sign
, j
+ sign
]
70 then ratio2 * magnitudesData.[i
+ sign
, j
+ sign
] + ratio1 * magnitudesData.[i
+ sign
, j
]
71 elif
angle < 3.f
* PI / 4.f
72 then ratio1 * magnitudesData.[i
+ sign
, j
] + ratio2 * magnitudesData.[i
+ sign
, j
- sign
]
74 then ratio2 * magnitudesData.[i
+ sign
, j
- sign
] + ratio1 * magnitudesData.[i
, j
- sign
]
75 elif
angle < 5.f
* PI / 4.f
76 then ratio1 * magnitudesData.[i
, j
- sign
] + ratio2 * magnitudesData.[i
- sign
, j
- sign
]
77 elif
angle < 3.f
* PI / 2.f
78 then ratio2 * magnitudesData.[i
- sign
, j
- sign
] + ratio1 * magnitudesData.[i
- sign
, j
]
79 elif
angle < 7.f
* PI / 4.f
80 then ratio1 * magnitudesData.[i
- sign
, j
] + ratio2 * magnitudesData.[i
- sign
, j
+ sign
]
81 else ratio2 * magnitudesData.[i
- sign
, j
+ sign
] + ratio1 * magnitudesData.[i
, j
+ sign
]
83 let m = magnitudesData.[i
, j
]
84 if m >= thresholdLow
&& m > mNeigbors 1 && m > mNeigbors -1
88 // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
90 IO.saveMat
magnitudes "magnitudes.png"
91 IO.saveMat
nms "nms.png"
93 let edges = new Matrix<byte
>(xGradient.Size)
94 let edgesData = edges.Data
96 // Hysteresis thresholding.
97 let toVisit = Stack<Point>()
98 for i
in 0 .. h - 1 do
99 for j
in 0 .. w - 1 do
100 if nmsData.[i
, j
] = 1uy && magnitudesData.[i
, j
] >= thresholdHigh
102 nmsData.[i
, j
] <- 0uy
103 toVisit.Push(Point(j
, i
))
104 while toVisit.Count > 0 do
105 let p = toVisit.Pop()
106 edgesData.[p.Y, p.X] <- 1uy
109 if i
' <> 0 || j' <> 0
113 if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
115 nmsData.[ni, nj] <- 0uy
116 toVisit.Push(Point(nj, ni))
118 edges, xGradient, yGradient