Parasitemia/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>(array2D [[ -1.0f; 0.0f; 1.0f ]
  29                                      [ -2.0f; 0.0f; 2.0f ]
  30                                      [ -1.0f; 0.0f; 1.0f ]])
  31
  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))
  36
  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.
  40
  41     let thresholdHigh, thresholdLow =
  42         let threshold, _, _ = otsu (histogramMat magnitudes 300)
  43         threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
  44
  45     // Non-maximum suppression.
  46     use nms = new Matrix<byte>(xGradient.Size)
  47
  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
  53
  54     for i = 1 to h - 2 do
  55         for j = 1 to w - 2 do
  56             let vx = xGradientData.[i, j]
  57             let vy = yGradientData.[i, j]
  58             if vx <> 0.f || vy <> 0.f
  59             then
  60                 let angle = anglesData.[i, j]
  61
  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
  65
  66                 let mNeigbors (sign: int) : float32 =
  67                     if angle < PI / 4.f
  68                     then ratio1 * magnitudesData.[i, j + sign] + ratio2 * magnitudesData.[i + sign, j + sign]
  69                     elif angle < PI / 2.f
  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]
  73                     elif angle < PI
  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]
  82
  83                 let m = magnitudesData.[i, j]
  84                 if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1
  85                 then
  86                     nmsData.[i, j] <- 1uy
  87
  88     // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
  89
  90     let edges = new Matrix<byte>(xGradient.Size)
  91     let edgesData = edges.Data
  92
  93     // Hysteresis thresholding.
  94     let toVisit = Stack<Point>()
  95     for i = 0 to h - 1 do
  96         for j = 0 to w - 1 do
  97             if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh
  98             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
 107                             then
 108                                 let ni = p.Y + i'
 109                                 let nj = p.X + j'
 110                                 if ni >= 0 && ni < h && nj >= 0 && nj < w && nmsData.[ni, nj] = 1uy
 111                                 then
 112                                     nmsData.[ni, nj] <- 0uy
 113                                     toVisit.Push(Point(nj, ni))
 114
 115     edges, xGradient, yGradient