1
module ParasitemiaCore.Edges
5 open System.Collections.Generic
15 /// Find edges of an image by using the Canny approach.
16 /// The thresholds are automatically defined with otsu on gradient magnitudes.
18 /// <param name="img"></param>
19 let find (img
: Image<Gray, float32
>) : Matrix<byte
> * Matrix<float32
> * Matrix<float32
> =
24 new Matrix<float32
>(array2D
[[ 1.0f; 0.0f; -1.0f ]
26 [ 1.0f; 0.0f; -1.0f ]])
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))
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.
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)
43 // Non-maximum suppression.
44 use nms = new Matrix<byte
>(xGradient.Size)
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
52 for i
in 0 .. h - 1 do
54 nmsData.[i
, w - 1] <- 0uy
56 for j
in 0 .. w - 1 do
58 nmsData.[h - 1, j
] <- 0uy
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
66 let angle = anglesData.[i
, j
]
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
72 let mNeigbors (sign
: int) : float32
=
74 then ratio1 * magnitudesData.[i
, j
+ sign
] + ratio2 * magnitudesData.[i
+ sign
, j
+ sign
]
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
]
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
]
89 let m = magnitudesData.[i
, j
]
90 if m >= thresholdLow
&& m > mNeigbors 1 && m > mNeigbors -1
94 // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
96 let edges = new Matrix<byte
>(xGradient.Size)
97 let edgesData = edges.Data
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
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
112 if i
' <> 0 || j' <> 0
116 if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
118 nmsData.[ni, nj] <- 0uy
119 toVisit.Push(Point(nj, ni))
121 edges, xGradient, yGradient