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>(
  29             array2D [[ -1.0f; 0.0f; 1.0f ]
  30                      [ -2.0f; 0.0f; 2.0f ]
  31                      [ -1.0f; 0.0f; 1.0f ]]
  32         )
  33
  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))
  38
  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.
  42
  43     let thresholdHigh, thresholdLow =
  44         let threshold, _, _ = otsu (histogramMat magnitudes 300)
  45         threshold + (sensibilityHigh * threshold), threshold - (sensibilityLow * threshold)
  46
  47     // Non-maximum suppression.
  48     use nms = new Matrix<byte>(xGradient.Size)
  49
  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
  55
  56     for i = 1 to h - 2 do
  57         for j = 1 to w - 2 do
  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]
  62
  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
  66
  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]
  74                     elif angle < PI then
  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]
  82                     else
  83                         ratio2 * magnitudesData.[i - sign, j + sign] + ratio1 * magnitudesData.[i, j + sign]
  84
  85                 let m = magnitudesData.[i, j]
  86                 if m >= thresholdLow && m > mNeigbors 1 && m > mNeigbors -1 then
  87                     nmsData.[i, j] <- 1uy
  88
  89     // suppressMConnections nms // It's not helpful for the rest of the process (ellipse detection).
  90
  91     let edges = new Matrix<byte>(xGradient.Size)
  92     let edgesData = edges.Data
  93
  94     // Hysteresis thresholding.
  95     let toVisit = Stack<Point>()
  96     for i = 0 to h - 1 do
  97         for j = 0 to w - 1 do
  98             if nmsData.[i, j] = 1uy && magnitudesData.[i, j] >= thresholdHigh 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 then
 107                                 let ni = p.Y + i'
 108                                 let nj = p.X + j'
 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))
 112
 113     edges, xGradient, yGradient